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

Changeset 3625


Ignore:
Timestamp:
2012-11-21T14:19:18+01:00 (10 years ago)
Author:
acc
Message:

Branch dev_NOC_2012_r3555. #1006. Step 7. Check in code now merged with dev_r3385_NOCS04_HAMF

Location:
branches/2012/dev_NOC_2012_rev3555
Files:
106 edited

Legend:

Unmodified
Added
Removed

  • branches/2012/dev_NOC_2012_rev3555/DOC/TexFiles/Chapters/Chap_SBC.tex

    r3609 r3625  
    11461146\label{SBC_cice} 
    11471147 
    1148 It is now possible to couple a global NEMO configuration (without AGRIF) to the CICE sea-ice 
     1148It is now possible to couple a regional or global NEMO configuration (without AGRIF) to the CICE sea-ice 
    11491149model by using \key{cice}.  The CICE code can be obtained from  
    11501150\href{http://oceans11.lanl.gov/trac/CICE/}{LANL} and the additional 'hadgem3' drivers will be required,  

  • branches/2012/dev_NOC_2012_rev3555/DOC/TexFiles/Chapters/Introduction.tex

    r3308 r3625  
    6363\citep{OASIS2006}. Two-way nesting is also available through an interface to the 
    6464AGRIF package (Adaptative Grid Refinement in \textsc{Fortran}) \citep{Debreu_al_CG2008}. 
    65 The interface code for coupling to an alternative sea ice model (CICE, \citet{Hunke2008}) is now  
    66 available although this is currently only designed for global domains, without the use of AGRIF. 
     65The interface code for coupling to an alternative sea ice model (CICE, \citet{Hunke2008}) 
     66has now been upgraded so that it works for both global and regional domains, although AGRIF  
     67is still not available. 
    6768 
    6869Other model characteristics are the lateral boundary conditions (chapter~\ref{LBC}).   

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/CONFIG/AMM12/EXP00/namelist

    r3609 r3625  
    137137                           !  =1 use observed ice-cover      , 
    138138                           !  =2 ice-model used                         ("key_lim3" or "key_lim2) 
     139   nn_ice_embd = 0         !  =0 levitating ice (no mass exchange, concentration/dilution effect) 
     140                           !  =1 levitating ice with mass and salt exchange but no presure effect 
     141                           !  =2 embedded sea-ice (full salt and mass exchanges and pressure) 
    139142   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    140143   ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/CONFIG/AMM12_PISCES/EXP00/namelist

    r3609 r3625  
    137137                           !  =1 use observed ice-cover      , 
    138138                           !  =2 ice-model used                         ("key_lim3" or "key_lim2) 
     139   nn_ice_embd = 0         !  =0 levitating ice (no mass exchange, concentration/dilution effect) 
     140                           !  =1 levitating ice with mass and salt exchange but no presure effect 
     141                           !  =2 embedded sea-ice (full salt and mass exchanges and pressure) 
    139142   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    140143   ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/CONFIG/GYRE/EXP00/namelist

    r3614 r3625  
    137137                           !  =1 use observed ice-cover      , 
    138138                           !  =2 ice-model used                         ("key_lim3" or "key_lim2) 
     139   nn_ice_embd = 0         !  =0 levitating ice (no mass exchange, concentration/dilution effect) 
     140                           !  =1 levitating ice with mass and salt exchange but no presure effect 
     141                           !  =2 embedded sea-ice (full salt and mass exchanges and pressure) 
    139142   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    140143   ln_rnf      = .false.   !  runoffs                                   (T => fill namsbc_rnf) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml

    r3294 r3625  
    124124 
    125125   <field id="empmr"        description="Net Upward Water Flux"                                        unit="kg/m2/s"  /> 
    126    <field id="empsmr"       description="concentration/dilution water flux"                            unit="kg/m2/s" /> 
     126   <field id="saltflx"      description="Downward Salt Flux"                                           unit="PSU/m2/s" /> 
    127127   <field id="snowpre"      description="Snow precipitation"                                           unit="kg/m2/s"  /> 
    128128   <field id="runoffs"      description="River Runoffs"                                                unit="Kg/m2/s"  /> 
    … …  
    145145   <field id="qsb_oce"      description="Sensible Downward Heat Flux over open ocean"                  unit="W/m2"     /> 
    146146   <field id="qla_oce"      description="Latent Downward Heat Flux over open ocean"                    unit="W/m2"     /> 
     147   <field id="qhc_oce"      description="Downward Heat Content of E-P over open ocean"                 unit="W/m2"     /> 
    147148   <field id="taum_oce"     description="wind stress module over open ocean"                           unit="N/m2"     /> 
    148149 
    … …  
    173174   <field id="v_imasstr"    description="Sea-ice mass transport along j-axis"                          unit="kg/s"     /> 
    174175 
     176   <!-- available if not defined key_vvl --> 
     177   <field id="emp_x_sst"      description="Concentration/Dilution term on SST"                         unit="kgC/m2/s" /> 
     178   <field id="emp_x_sss"      description="Concentration/Dilution term on SSS"                       unit="kgPSU/m2/s" /> 
    175179   <!-- available key_coupled --> 
    176180   <field id="snow_ao_cea"  description="Snow over ice-free ocean (cell average)"                      unit="kg/m2/s"  /> 
    … …  
    10161020     <field ref="empmr"        name="sowaflup"  /> 
    10171021     <field ref="qsr"          name="soshfldo"  /> 
    1018      <field ref="empsmr"       name="sowaflcd"  /> 
     1022     <field ref="saltflx"      name="sosfldow"  /> 
    10191023     <field ref="qt"           name="sohefldo"  /> 
    10201024     <field ref="mldr10_1"     name="somxl010"  /> 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist

    r3609 r3625  
    137137                           !  =1 use observed ice-cover      , 
    138138                           !  =2 ice-model used                         ("key_lim3" or "key_lim2) 
     139   nn_ice_embd = 0         !  =0 levitating ice (no mass exchange, concentration/dilution effect) 
     140                           !  =1 levitating ice with mass and salt exchange but no presure effect 
     141                           !  =2 embedded sea-ice (full salt and mass exchanges and pressure) 
    139142   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    140143   ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist

    r3614 r3625  
    137137                           !  =1 use observed ice-cover      , 
    138138                           !  =2 ice-model used                         ("key_lim3" or "key_lim2) 
     139   nn_ice_embd = 0         !  =0 levitating ice (no mass exchange, concentration/dilution effect) 
     140                           !  =1 levitating ice with mass and salt exchange but no presure effect 
     141                           !  =2 embedded sea-ice (full salt and mass exchanges and pressure) 
    139142   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    140143   ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf) 
    … …  
    650653   sn_mld  = 'dyna_grid_T' ,    120            , 'somixhgt' ,  .true.    , .true. ,   'yearly'  , ''       , '' 
    651654   sn_emp  = 'dyna_grid_T' ,    120            , 'sowaflcd' ,  .true.    , .true. ,   'yearly'  , ''       , '' 
     655!  sn_emp  = 'dyna_grid_T' ,    120            , 'sowaflup' ,  .true.    , .true. ,   'yearly'  , ''       , '' ! v3.5+ 
     656!  sn_sfx  = 'dyna_grid_T' ,    120            , 'sosfldow' ,  .true.    , .true. ,   'yearly'  , ''       , '' ! v3.5+ 
    652657   sn_ice  = 'dyna_grid_T' ,    120            , 'soicecov' ,  .true.    , .true. ,   'yearly'  , ''       , '' 
    653658   sn_qsr  = 'dyna_grid_T' ,    120            , 'soshfldo' ,  .true.    , .true. ,   'yearly'  , ''       , '' 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90

    r2715 r3625  
    1919   PUBLIC    ice_alloc_2  !  Called in iceini_2.F90 
    2020 
    21    INTEGER , PUBLIC ::   numit     !: ice iteration index 
    22    REAL(wp), PUBLIC ::   rdt_ice   !: ice time step 
     21   INTEGER , PUBLIC ::   numit        !: ice iteration index 
     22   REAL(wp), PUBLIC ::   rdt_ice      !: ice time step 
    2323 
    2424   !                                                                     !!* namicerun read in iceini  * 
    … …  
    9898   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qstoif        !: Energy stored in the brine pockets 
    9999   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fbif          !: Heat flux at the ice base 
    100    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmsnif       !: Variation of snow mass 
    101    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmicif       !: Variation of ice mass 
     100   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_snw       !: Variation of snow mass over 1 time step           [Kg/m2] 
     101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_snw       !: Heat content associated with rdm_snw              [J/m2] 
     102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_ice       !: Variation of ice  mass over 1 time step           [Kg/m2] 
     103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_ice       !: Heat content associated with rdm_ice              [J/m2] 
    102104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qldif         !: heat balance of the lead (or of the open ocean) 
    103105   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qcmif         !: Energy needed to freeze the ocean surface layer 
    … …  
    153155 
    154156      ALLOCATE(phicif(jpi,jpj) , pfrld  (jpi,jpj) , qstoif (jpi,jpj) ,     & 
    155          &     fbif  (jpi,jpj) , rdmsnif(jpi,jpj) , rdmicif(jpi,jpj) ,     & 
     157         &     fbif  (jpi,jpj) , rdm_snw(jpi,jpj) , rdq_snw(jpi,jpj) ,     & 
     158         &                       rdm_ice(jpi,jpj) , rdq_ice(jpi,jpj) ,     & 
    156159         &     qldif (jpi,jpj) , qcmif  (jpi,jpj) , fdtcn  (jpi,jpj) ,     & 
    157160         &     qdtcn (jpi,jpj) , thcm   (jpi,jpj)                    , STAT=ierr(4) ) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90

    r3294 r3625  
    1313   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    15    !!   ice_init_2       : sea-ice model initialization 
    16    !!   ice_run_2        : Definition some run parameter for ice model 
     15   !!   ice_init_2    : sea-ice model initialization 
     16   !!   ice_run_2     : Definition some run parameter for ice model 
    1717   !!---------------------------------------------------------------------- 
    18    USE phycst           ! physical constants 
    19    USE dom_oce          ! ocean domain 
    20    USE sbc_oce          ! surface boundary condition: ocean 
    21    USE sbc_ice          ! LIM2 surface boundary condition 
    22    USE dom_ice_2        ! LIM2 ice domain 
    23    USE par_ice_2        ! LIM2 parameters 
    24    USE thd_ice_2        ! LIM2 thermodynamical variables 
    25    USE ice_2            ! LIM2 ice variable 
    26    USE limmsh_2         ! LIM2 mesh 
    27    USE limistate_2      ! LIM2 initial state 
    28    USE limrst_2         ! LIM2 restart 
    29    USE limsbc_2         ! LIM2 surface boundary condition 
    30    USE in_out_manager   ! I/O manager 
    31    USE lib_mpp          ! MPP library 
     18   USE phycst         ! physical constants 
     19   USE dom_oce        ! ocean domain 
     20   USE sbc_oce        ! surface boundary condition: ocean 
     21   USE sbc_ice        ! LIM2 surface boundary condition 
     22   USE dom_ice_2      ! LIM2 ice domain 
     23   USE par_ice_2      ! LIM2 parameters 
     24   USE thd_ice_2      ! LIM2 thermodynamical variables 
     25   USE ice_2          ! LIM2 ice variable 
     26   USE limmsh_2       ! LIM2 mesh 
     27   USE limistate_2    ! LIM2 initial state 
     28   USE limrst_2       ! LIM2 restart 
     29   USE limsbc_2       ! LIM2 surface boundary condition 
     30   USE in_out_manager ! I/O manager 
     31   USE lib_mpp        ! MPP library 
     32   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3233 
    3334   IMPLICIT NONE 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90

    r3294 r3625  
    1414   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    1515   !!---------------------------------------------------------------------- 
    16    !!   lim_adv_x_2  : advection of sea ice on x axis 
    17    !!   lim_adv_y_2  : advection of sea ice on y axis 
     16   !!   lim_adv_x_2   : advection of sea ice on x axis 
     17   !!   lim_adv_y_2   : advection of sea ice on y axis 
    1818   !!---------------------------------------------------------------------- 
    1919   USE dom_oce 
    … …  
    2121   USE ice_2 
    2222   USE lbclnk 
    23    USE in_out_manager     ! I/O manager 
    24    USE lib_mpp            ! MPP library 
    25    USE wrk_nemo           ! work arrays 
    26    USE prtctl             ! Print control 
     23   USE in_out_manager ! I/O manager 
     24   USE lib_mpp        ! MPP library 
     25   USE wrk_nemo       ! work arrays 
     26   USE prtctl         ! Print control 
     27   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2728 
    2829   IMPLICIT NONE 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdia_2.F90

    r2715 r3625  
    2424   USE in_out_manager  ! I/O manager 
    2525   USE lib_mpp         ! MPP library 
     26   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2627 
    2728   IMPLICIT NONE 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90

    r2715 r3625  
    1919   USE in_out_manager  ! I/O manager 
    2020   USE lib_mpp         ! MPP library 
     21   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2122 
    2223   IMPLICIT NONE 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90

    r3294 r3625  
    3131   USE in_out_manager   ! I/O manager 
    3232   USE prtctl           ! Print control 
     33   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3334 
    3435   IMPLICIT NONE 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90

    r3294 r3625  
    2121   USE prtctl           ! Print control 
    2222   USE in_out_manager   ! I/O manager 
     23   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2324 
    2425   IMPLICIT NONE 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90

    r3294 r3625  
    2727   USE iom 
    2828   USE in_out_manager 
     29   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2930 
    3031   IMPLICIT NONE 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90

    r3294 r3625  
    2323   USE wrk_nemo         ! work arrays 
    2424#endif 
     25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2526 
    2627   IMPLICIT NONE 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90

    r3294 r3625  
    3030   USE in_out_manager ! I/O manager 
    3131   USE prtctl         ! Print control 
     32   USE oce     , ONLY : snwice_mass, snwice_mass_b 
     33   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3234 
    3335   IMPLICIT NONE 
    … …  
    8082      REAL(wp) ::   zs21_11, zs21_12, zs21_21, zs21_22 
    8183      REAL(wp) ::   zs22_11, zs22_12, zs22_21, zs22_22 
     84      REAL(wp) ::   zintb, zintn 
    8285      REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld, zmass, zcorl 
    8386      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct, za2ct, zresr 
    8487      REAL(wp), POINTER, DIMENSION(:,:) ::   zc1u, zc1v, zc2u, zc2v 
    85       REAL(wp), POINTER, DIMENSION(:,:) ::   zsang 
     88      REAL(wp), POINTER, DIMENSION(:,:) ::   zsang, zpice 
    8689      REAL(wp), POINTER, DIMENSION(:,:) ::   zu0, zv0 
    8790      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_n, zv_n 
    … …  
    9396       
    9497      CALL wrk_alloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr ) 
    95       CALL wrk_alloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang ) 
     98      CALL wrk_alloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice ) 
    9699      CALL wrk_alloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 ) 
    97100      CALL wrk_alloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 ) 
    … …  
    129132!i    zviszeta(:,jpj+1) = 0._wp    ;    zviseta(:,jpj+1) = 0._wp 
    130133 
     134      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
     135          ! 
     136          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
     137          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 
     138         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 
     139          ! 
     140          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
     141          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
     142         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
     143          ! 
     144         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0 
     145          ! 
     146         ! 
     147      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
     148         zpice(:,:) = ssh_m(:,:) 
     149      ENDIF 
    131150 
    132151      ! Ice mass, ice strength, and wind stress at the center            | 
    … …  
    196215 
    197216            ! Gradient of the sea surface height 
    198             zgsshx =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji-1,jj  ))/e1u(ji-1,jj  )   & 
    199                &       +  (ssh_m(ji  ,jj-1) - ssh_m(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp 
    200             zgsshy =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji  ,jj-1))/e2v(ji  ,jj-1)   & 
    201                &       +  (ssh_m(ji-1,jj  ) - ssh_m(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp 
     217            zgsshx =  (   (zpice(ji  ,jj  ) - zpice(ji-1,jj  ))/e1u(ji-1,jj  )   & 
     218               &       +  (zpice(ji  ,jj-1) - zpice(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp 
     219            zgsshy =  (   (zpice(ji  ,jj  ) - zpice(ji  ,jj-1))/e2v(ji  ,jj-1)   & 
     220               &       +  (zpice(ji-1,jj  ) - zpice(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp 
    202221 
    203222            ! Computation of the velocity field taking into account the ice-ice interaction.                                  
    … …  
    575594 
    576595      CALL wrk_dealloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr ) 
    577       CALL wrk_dealloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang ) 
     596      CALL wrk_dealloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice ) 
    578597      CALL wrk_dealloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 ) 
    579598      CALL wrk_dealloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 ) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90

    r3294 r3625  
    99   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case 
    1010   !!             -   ! 2010-11 (G. Madec) ice-ocean stress computed at each ocean time-step 
    11    !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation 
     11   !!           3.3.1 ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation 
     12   !!            3.5  ! 2012-11 ((G. Madec, Y. Aksenov, A. Coward) salt and heat fluxes associated with e-p 
    1213   !!---------------------------------------------------------------------- 
    1314#if defined key_lim2 
    … …  
    2829   USE sbc_oce          ! surface boundary condition: ocean 
    2930   USE sbccpl 
    30  
     31   USE cpl_oasis3, ONLY : lk_cpl 
     32   USE oce       , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass  
    3133   USE albedo           ! albedo parameters 
    3234   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges 
    … …  
    3739   USE iom              ! I/O library 
    3840   USE prtctl           ! Print control 
    39    USE cpl_oasis3, ONLY : lk_cpl 
     41   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    4042 
    4143   IMPLICIT NONE 
    … …  
    8890      !!              - Update the fluxes provided to the ocean 
    8991      !!      
    90       !! ** Outputs : - qsr     : sea heat flux:    solar  
    91       !!              - qns     : sea heat flux: non solar 
    92       !!              - emp     : freshwater budget: volume flux  
    93       !!              - emps    : freshwater budget: concentration/dillution  
     92      !! ** Outputs : - qsr     : sea heat flux    : solar  
     93      !!              - qns     : sea heat flux    : non solar (including heat content of the mass flux) 
     94      !!              - emp     : freshwater budget: mass flux  
     95      !!              - sfx     : freshwater budget: salt flux due to Freezing/Melting 
    9496      !!              - utau    : sea surface i-stress (ocean referential) 
    9597      !!              - vtau    : sea surface j-stress (ocean referential) 
    … …  
    107109      INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       - 
    108110      INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       - 
    109       REAL(wp) ::   zqsr, zqns, zfm            ! local scalars 
    110       REAL(wp) ::   zinda, zfons, zemp         !   -      - 
     111      REAL(wp) ::   zqsr,     zqns,   zfmm     ! local scalars 
     112      REAL(wp) ::   zinda,    zfsalt, zemp     !   -      - 
     113      REAL(wp) ::   zemp_snw, zqhc,   zcd      !   -      - 
     114      REAL(wp) ::   zswitch                    !   -      - 
    111115      REAL(wp), POINTER, DIMENSION(:,:)   ::   zqnsoce       ! 2D workspace 
    112116      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace 
    … …  
    115119      CALL wrk_alloc( jpi, jpj, zqnsoce ) 
    116120      CALL wrk_alloc( jpi, jpj, 1, zalb, zalbp ) 
     121 
     122      SELECT CASE( nn_ice_embd )                 ! levitating or embedded sea-ice option 
     123        CASE( 0    )   ;   zswitch = 1           ! (0) standard levitating sea-ice : salt exchange only 
     124        CASE( 1, 2 )   ;   zswitch = 0           ! (1) levitating sea-ice: salt and volume exchange but no pressure effect 
     125                                                 ! (2) embedded sea-ice : salt and volume fluxes and pressure 
     126      END SELECT                                 !     
    117127 
    118128      !------------------------------------------! 
    … …  
    133143            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv  
    134144 
    135 !!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   !   = 0. if pure ocean else 1. (at previous time) 
    136 !!$ 
    137 !!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   !   = 0. if pure ocean else 1. (at current  time) 
    138 !!$ 
    139 !!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      !   = 1. if (snow and no ice at previous time) else 0. ??? 
    140 !!$            ELSE                             ;   ifvt = 0. 
     145!!$            attempt to explain the tricky flags set above.... 
     146!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice adv, but before ice thermo) 
     147!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice thermo) 
     148!!$ 
     149!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      ! = zinda if previous thermodynamic step overmelted the ice??? 
     150!!$            ELSE                             ;   ifvt = 0.         !  
    141151!!$            ENDIF 
    142152!!$ 
    143 !!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases from previous to current 
     153!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases due to ice thermodynamics 
    144154!!$            ELSE                                     ;   idfr = 1.    
    145155!!$            ENDIF 
    146156!!$ 
    147 !!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous and pure ocean at current 
     157!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous time and ice-free ocean currently 
    148158!!$ 
    149159!!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr 
     160!!$                    = i1mfr if ifvt = 1 i.e.  
     161!!$                    = idfr  if ifvt = 0 
    150162!!$!                 snow no ice   ice         ice or nothing  lead fraction increases 
    151163!!$!                 at previous   now           at previous 
    152 !!$!                -> ice aera increases  ???         -> ice aera decreases ??? 
     164!!$!                -> ice area increases  ???         -> ice area decreases ??? 
    153165!!$ 
    154166!!$            iadv    = ( 1  - i1mfr ) * zinda 
    … …  
    174186#endif             
    175187            !  computation the non solar heat flux at ocean surface 
    176             zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads 
    177                &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            & 
    178                &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice    & 
    179                &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice  
    180  
    181             fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ??? 
    182             ! 
     188            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr                                              &   ! part of the solar energy used in leads 
     189               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                             & 
     190               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice  & 
     191               &       + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice  
     192 
     193            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! store residual heat flux (to put into the ocean at the next time-step) 
     194            zqhc = ( rdq_snw(ji,jj)                                     & 
     195                 & + rdq_ice(ji,jj) * ( 1.- zswitch) ) * r1_rdtice       ! heat flux due to snow ( & ice heat content,  
     196            !                                                            !           if ice/ocean mass exchange active)  
    183197            qsr  (ji,jj) = zqsr                                          ! solar heat flux  
    184             qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux 
     198            qns  (ji,jj) = zqns - fdtcn(ji,jj) + zqhc                    ! non solar heat flux  
     199            ! 
     200            !                          !------------------------------------------! 
     201            !                          !  mass and salt flux at the ocean surface ! 
     202            !                          !------------------------------------------! 
     203            ! 
     204            ! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area) 
     205#if defined key_coupled 
     206            !                                                  ! coupled mode:  
     207            zemp = + emp_tot(ji,jj)                            &     ! net mass flux over the grid cell (ice+ocean area) 
     208               &   - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )          ! minus the mass flux intercepted by sea-ice 
     209#else 
     210            !                                                  ! forced  mode:  
     211            zemp = + emp(ji,jj)     *         frld(ji,jj)      &     ! mass flux over open ocean fraction  
     212               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &     ! liquid precip. over ice reaches directly the ocean 
     213               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )          ! snow is intercepted by sea-ice (previous frld) 
     214#endif             
     215            ! 
     216            ! mass flux at the ocean/ice interface (sea ice fraction) 
     217            zemp_snw = rdm_snw(ji,jj) * r1_rdtice                    ! snow melting = pure water that enters the ocean 
     218            zfmm     = rdm_ice(ji,jj) * r1_rdtice                    ! Freezing minus Melting (F-M) 
     219 
     220            ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s] 
     221            zfsalt = - sice_0(ji,jj) * zfmm                          ! F-M salt exchange 
     222            zcd    =   soce_0(ji,jj) * zfmm                          ! concentration/dilution term due to F-M 
     223            ! 
     224            ! salt flux only       : add concentration dilution term in salt flux  and no  F-M term in volume flux 
     225            ! salt and mass fluxes : non concentration dilution term in salt flux  and add F-M term in volume flux 
     226            sfx (ji,jj) = zfsalt +                  zswitch  * zcd   ! salt flux (+ C/D if no ice/ocean mass exchange) 
     227            emp (ji,jj) = zemp   + zemp_snw + ( 1.- zswitch) * zfmm  ! mass flux (+ F/M mass flux if ice/ocean mass exchange) 
     228            ! 
    185229         END DO 
    186230      END DO 
     231      !                                !------------------------------------------! 
     232      !                                !    mass of snow and ice per unit area    ! 
     233      !                                !------------------------------------------! 
     234      IF( nn_ice_embd /= 0 ) THEN      ! embedded sea-ice (mass required) 
     235         snwice_mass_b(:,:) = snwice_mass(:,:)                  ! save mass from the previous ice time step 
     236         !                                                      ! new mass per unit area 
     237         snwice_mass  (:,:) = tms(:,:) * ( rhosn * hsnif(:,:) + rhoic * hicif(:,:)  ) * ( 1.0 - frld(:,:) ) 
     238         !                                                      ! time evolution of snow+ice mass 
     239         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / rdt_ice 
     240      ENDIF 
    187241 
    188242      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )       
    … …  
    190244      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1.e0 - pfrld(:,:)) ) 
    191245 
    192       !------------------------------------------! 
    193       !      mass flux at the ocean surface      ! 
    194       !------------------------------------------! 
    195       DO jj = 1, jpj 
    196          DO ji = 1, jpi 
    197             ! 
    198 #if defined key_coupled 
    199             ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode) 
    200             zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  
    201                &   + rdmsnif(ji,jj) * r1_rdtice                                   !  freshwaterflux due to snow melting  
    202 #else 
    203             !  computing freshwater exchanges at the ice/ocean interface 
    204             zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction  
    205                &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean 
    206                &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  change in ice cover within the time step 
    207                &   + rdmsnif(ji,jj) * r1_rdtice                    !  freshwater flux due to snow melting  
    208 #endif             
    209             ! 
    210             !  computing salt exchanges at the ice/ocean interface 
    211             zfons = ( soce_0(ji,jj) - sice_0(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice )  
    212             ! 
    213             !  converting the salt flux from ice to a freshwater flux from ocean 
    214             zfm  = zfons / ( sss_m(ji,jj) + epsi16 ) 
    215             ! 
    216             emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution) 
    217             emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution) 
    218             ! 
    219          END DO 
    220       END DO 
    221  
    222246      IF( lk_diaar5 ) THEN       ! AR5 diagnostics 
    223          CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * r1_rdtice ) 
    224          CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdmicif(:,:) * r1_rdtice ) 
    225          CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdmicif(:,:) * r1_rdtice ) 
     247         CALL iom_put( 'isnwmlt_cea'  ,                 rdm_snw(:,:) * r1_rdtice ) 
     248         CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdm_ice(:,:) * r1_rdtice ) 
     249         CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdm_ice(:,:) * r1_rdtice ) 
    226250      ENDIF 
    227251 
    … …  
    243267      IF(ln_ctl) THEN            ! control print 
    244268         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ') 
    245          CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ') 
     269         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=sfx   , clinfo2=' sfx     : ') 
    246270         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   & 
    247271            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask ) 
    … …  
    439463         END WHERE 
    440464      ENDIF 
     465      !                                      ! embedded sea ice 
     466      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass 
     467         snwice_mass  (:,:) = tms(:,:) * ( rhosn * hsnif(:,:) + rhoic * hicif(:,:)  ) * ( 1.0 - frld(:,:) ) 
     468         snwice_mass_b(:,:) = snwice_mass(:,:) 
     469      ELSE 
     470         snwice_mass  (:,:) = 0.e0           ! no mass exchanges 
     471         snwice_mass_b(:,:) = 0.e0           ! no mass exchanges 
     472      ENDIF 
     473      IF( nn_ice_embd == 2 .AND.          &  ! full embedment (case 2) & no restart :  
     474         &   .NOT.ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area 
     475         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
     476         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 
     477      ENDIF 
    441478      ! 
    442479   END SUBROUTINE lim_sbc_init_2 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90

    r3294 r3625  
    1313   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    15    !!   lim_thd_2      : thermodynamic of sea ice 
    16    !!   lim_thd_init_2 : initialisation of sea-ice thermodynamic 
     15   !!   lim_thd_2       : thermodynamic of sea ice 
     16   !!   lim_thd_init_2  : initialisation of sea-ice thermodynamic 
    1717   !!---------------------------------------------------------------------- 
    18    USE phycst          ! physical constants 
    19    USE dom_oce         ! ocean space and time domain variables 
     18   USE phycst           ! physical constants 
     19   USE dom_oce          ! ocean space and time domain variables 
    2020   USE domvvl 
    2121   USE lbclnk 
    22    USE in_out_manager  ! I/O manager 
     22   USE in_out_manager   ! I/O manager 
    2323   USE lib_mpp 
    24    USE wrk_nemo        ! work arrays 
    25    USE iom             ! IOM library 
    26    USE ice_2           ! LIM sea-ice variables 
    27    USE sbc_oce         !  
    28    USE sbc_ice         !  
    29    USE thd_ice_2       ! LIM thermodynamic sea-ice variables 
    30    USE dom_ice_2       ! LIM sea-ice domain 
     24   USE wrk_nemo         ! work arrays 
     25   USE iom              ! IOM library 
     26   USE ice_2            ! LIM sea-ice variables 
     27   USE sbc_oce          !  
     28   USE sbc_ice          !  
     29   USE thd_ice_2        ! LIM thermodynamic sea-ice variables 
     30   USE dom_ice_2        ! LIM sea-ice domain 
    3131   USE limthd_zdf_2 
    3232   USE limthd_lac_2 
    3333   USE limtab_2 
    34    USE prtctl          ! Print control 
    35    USE cpl_oasis3, ONLY : lk_cpl 
    36    USE diaar5, ONLY :   lk_diaar5 
    37        
     34   USE prtctl           ! Print control 
     35   USE cpl_oasis3, ONLY :   lk_cpl 
     36   USE diaar5    , ONLY :   lk_diaar5 
     37   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     38    
    3839   IMPLICIT NONE 
    3940   PRIVATE 
    … …  
    5556   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    5657   !!---------------------------------------------------------------------- 
    57  
    5858CONTAINS 
    5959 
    … …  
    8989      REAL(wp) ::   za , zh, zthsnice    ! 
    9090      REAL(wp) ::   zfric_u              ! friction velocity  
    91       REAL(wp) ::   zfnsol               ! total non solar heat 
    92       REAL(wp) ::   zfontn               ! heat flux from snow thickness 
    9391      REAL(wp) ::   zfntlat, zpareff     ! test. the val. of lead heat budget 
    9492 
    … …  
    129127      zdvolif(:,:) = 0.e0   ! total variation of ice volume 
    130128      zdvonif(:,:) = 0.e0   ! transformation of snow to sea-ice volume 
    131 !      zdvonif(:,:) = 0.e0   ! lateral variation of ice volume 
    132129      zlicegr(:,:) = 0.e0   ! lateral variation of ice volume 
    133130      zdvomif(:,:) = 0.e0   ! variation of ice volume at bottom due to melting only 
    … …  
    137134      ffltbif(:,:) = 0.e0   ! linked with fstric 
    138135      qfvbq  (:,:) = 0.e0   ! linked with fstric 
    139       rdmsnif(:,:) = 0.e0   ! variation of snow mass per unit area 
    140       rdmicif(:,:) = 0.e0   ! variation of ice mass per unit area 
     136      rdm_snw(:,:) = 0.e0   ! variation of snow mass over 1 time step 
     137      rdq_snw(:,:) = 0.e0   ! heat content associated with rdm_snw 
     138      rdm_ice(:,:) = 0.e0   ! variation of ice mass over 1 time step 
     139      rdq_ice(:,:) = 0.e0   ! heat content associated with rdm_ice 
    141140      zmsk (:,:,:) = 0.e0 
    142141 
    … …  
    199198      !-------------------------------------------------------------------------- 
    200199 
    201       sst_m(:,:) = sst_m(:,:) + rt0 
    202  
    203 !CDIR NOVERRCHK 
    204       DO jj = 1, jpj 
    205 !CDIR NOVERRCHK 
     200      !CDIR NOVERRCHK 
     201      DO jj = 1, jpj 
     202         !CDIR NOVERRCHK 
    206203         DO ji = 1, jpi 
    207204            zthsnice       = hsnif(ji,jj) + hicif(ji,jj) 
    … …  
    217214            !  temperature and turbulent mixing (McPhee, 1992) 
    218215            zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )  ! friction velocity 
    219             fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) )  
     216            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) + rt0 - tfu(ji,jj) )  
    220217            qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice 
    221218                         
    222219            !  partial computation of the lead energy budget (qldif) 
    223220#if defined key_coupled  
    224             qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                             & 
     221            qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                                  & 
    225222               &    * (   ( qsr_tot(ji,jj) - qsr_ice(ji,jj,1) * zfricp ) * ( 1.0 - thcm(ji,jj) )   & 
    226223               &        + ( qns_tot(ji,jj) - qns_ice(ji,jj,1) * zfricp )                           & 
    227224               &        + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) )   ) 
    228225#else 
    229             zfontn         = ( sprecip(ji,jj) / rhosn ) * xlsn  !   energy for melting solid precipitation 
    230             zfnsol         = qns(ji,jj)                         !  total non solar flux over the ocean 
    231             qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )   & 
    232                &                               + zfnsol + fdtcn(ji,jj) - zfontn     & 
    233                &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   & 
    234                &                        * frld(ji,jj) * rdt_ice     
    235 !!$            qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj)  
    236 !!$               &           * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )      & 
    237 !!$               &             + qns(ji,jj)  + fdtcn(ji,jj) - zfontn     & 
    238 !!$               &             + ( 1.0 - zindb ) * fsbbq(ji,jj)      )   & 
     226            qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj)                    & 
     227               &                        * (  qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )    & 
     228               &                           + qns(ji,jj)  +  fdtcn(ji,jj)           & 
     229               &                           + ( 1.0 - zindb ) * fsbbq(ji,jj)      ) 
    239230#endif 
    240231            !  parlat : percentage of energy used for lateral ablation (0.0)  
    … …  
    246237             
    247238            !  energy needed to bring ocean surface layer until its freezing 
    248             qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1)   & 
    249                 &          * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda ) 
     239            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda ) 
    250240             
    251241            !  calculate oceanic heat flux. 
    … …  
    257247      END DO 
    258248       
    259       sst_m(:,:) = sst_m(:,:) - rt0 
    260                 
    261249      !         Select icy points and fulfill arrays for the vectorial grid. 
    262250      !---------------------------------------------------------------------- 
    … …  
    312300         CALL tab_2d_1d_2( nbpb, qldif_1d   (1:nbpb)     , qldif      , jpi, jpj, npb(1:nbpb) ) 
    313301         CALL tab_2d_1d_2( nbpb, qstbif_1d  (1:nbpb)     , qstoif     , jpi, jpj, npb(1:nbpb) ) 
    314          CALL tab_2d_1d_2( nbpb, rdmicif_1d (1:nbpb)     , rdmicif    , jpi, jpj, npb(1:nbpb) ) 
     302         CALL tab_2d_1d_2( nbpb, rdm_ice_1d (1:nbpb)     , rdm_ice    , jpi, jpj, npb(1:nbpb) ) 
     303         CALL tab_2d_1d_2( nbpb, rdq_ice_1d (1:nbpb)     , rdq_ice    , jpi, jpj, npb(1:nbpb) ) 
    315304         CALL tab_2d_1d_2( nbpb, dmgwi_1d   (1:nbpb)     , dmgwi      , jpi, jpj, npb(1:nbpb) ) 
     305         CALL tab_2d_1d_2( nbpb, rdm_snw_1d (1:nbpb)     , rdm_snw    , jpi, jpj, npb(1:nbpb) ) 
     306         CALL tab_2d_1d_2( nbpb, rdq_snw_1d (1:nbpb)     , rdq_snw    , jpi, jpj, npb(1:nbpb) ) 
    316307         CALL tab_2d_1d_2( nbpb, qlbbq_1d   (1:nbpb)     , zqlbsbq    , jpi, jpj, npb(1:nbpb) ) 
    317308         ! 
    … …  
    332323         CALL tab_1d_2d_2( nbpb, qfvbq      , npb, qfvbq_1d  (1:nbpb)     , jpi, jpj ) 
    333324         CALL tab_1d_2d_2( nbpb, qstoif     , npb, qstbif_1d (1:nbpb)     , jpi, jpj ) 
    334          CALL tab_1d_2d_2( nbpb, rdmicif    , npb, rdmicif_1d(1:nbpb)     , jpi, jpj ) 
     325         CALL tab_1d_2d_2( nbpb, rdm_ice    , npb, rdm_ice_1d(1:nbpb)     , jpi, jpj ) 
     326         CALL tab_1d_2d_2( nbpb, rdq_ice    , npb, rdq_ice_1d(1:nbpb)     , jpi, jpj ) 
    335327         CALL tab_1d_2d_2( nbpb, dmgwi      , npb, dmgwi_1d  (1:nbpb)     , jpi, jpj ) 
    336          CALL tab_1d_2d_2( nbpb, rdmsnif    , npb, rdmsnif_1d(1:nbpb)     , jpi, jpj ) 
     328         CALL tab_1d_2d_2( nbpb, rdm_snw    , npb, rdm_snw_1d(1:nbpb)     , jpi, jpj ) 
     329         CALL tab_1d_2d_2( nbpb, rdq_snw    , npb, rdq_snw_1d(1:nbpb)     , jpi, jpj ) 
    337330         CALL tab_1d_2d_2( nbpb, zdvosif    , npb, dvsbq_1d  (1:nbpb)     , jpi, jpj ) 
    338331         CALL tab_1d_2d_2( nbpb, zdvobif    , npb, dvbbq_1d  (1:nbpb)     , jpi, jpj ) 
    … …  
    393386      IF( nbpac > 0 ) THEN 
    394387         ! 
    395          zlicegr(:,:) = rdmicif(:,:)      ! to output the lateral sea-ice growth  
     388         zlicegr(:,:) = rdm_ice(:,:)      ! to output the lateral sea-ice growth  
    396389         !...Put the variable in a 1-D array for lateral accretion 
    397390         CALL tab_2d_1d_2( nbpac, frld_1d   (1:nbpac)     , frld       , jpi, jpj, npac(1:nbpac) ) 
    … …  
    404397         CALL tab_2d_1d_2( nbpac, qcmif_1d  (1:nbpac)     , qcmif      , jpi, jpj, npac(1:nbpac) ) 
    405398         CALL tab_2d_1d_2( nbpac, qstbif_1d (1:nbpac)     , qstoif     , jpi, jpj, npac(1:nbpac) ) 
    406          CALL tab_2d_1d_2( nbpac, rdmicif_1d(1:nbpac)     , rdmicif    , jpi, jpj, npac(1:nbpac) ) 
     399         CALL tab_2d_1d_2( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice    , jpi, jpj, npac(1:nbpac) ) 
     400         CALL tab_2d_1d_2( nbpac, rdq_ice_1d(1:nbpac)     , rdq_ice    , jpi, jpj, npac(1:nbpac) ) 
    407401         CALL tab_2d_1d_2( nbpac, dvlbq_1d  (1:nbpac)     , zdvolif    , jpi, jpj, npac(1:nbpac) ) 
    408402         CALL tab_2d_1d_2( nbpac, tfu_1d    (1:nbpac)     , tfu        , jpi, jpj, npac(1:nbpac) ) 
    … …  
    418412         CALL tab_1d_2d_2( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d   (1:nbpac , 3 ), jpi, jpj ) 
    419413         CALL tab_1d_2d_2( nbpac, qstoif     , npac(1:nbpac), qstbif_1d (1:nbpac)     , jpi, jpj ) 
    420          CALL tab_1d_2d_2( nbpac, rdmicif    , npac(1:nbpac), rdmicif_1d(1:nbpac)     , jpi, jpj ) 
     414         CALL tab_1d_2d_2( nbpac, rdm_ice    , npac(1:nbpac), rdm_ice_1d(1:nbpac)     , jpi, jpj ) 
     415         CALL tab_1d_2d_2( nbpac, rdq_ice    , npac(1:nbpac), rdq_ice_1d(1:nbpac)     , jpi, jpj ) 
    421416         CALL tab_1d_2d_2( nbpac, zdvolif    , npac(1:nbpac), dvlbq_1d  (1:nbpac)     , jpi, jpj ) 
    422417         ! 
    … …  
    449444      CALL iom_put( 'iceprod_cea' , hicifp (:,:) * zztmp     )   ! Ice produced               [m/s] 
    450445      IF( lk_diaar5 ) THEN 
    451          CALL iom_put( 'snowmel_cea' , rdmsnif(:,:) * zztmp     )   ! Snow melt                  [kg/m2/s] 
     446         CALL iom_put( 'snowmel_cea' , rdm_snw(:,:) * zztmp     )   ! Snow melt                  [kg/m2/s] 
    452447         zztmp = rhoic / rdt_ice 
    453448         CALL iom_put( 'sntoice_cea' , zdvonif(:,:) * zztmp     )   ! Snow to Ice transformation [kg/m2/s] 
    454449         CALL iom_put( 'ticemel_cea' , zdvosif(:,:) * zztmp     )   ! Melt at Sea Ice top        [kg/m2/s] 
    455450         CALL iom_put( 'bicemel_cea' , zdvomif(:,:) * zztmp     )   ! Melt at Sea Ice bottom     [kg/m2/s] 
    456          zlicegr(:,:) = MAX( 0.e0, rdmicif(:,:)-zlicegr(:,:) ) 
    457          CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Latereal sea ice growth    [kg/m2/s] 
     451         zlicegr(:,:) = MAX( 0.e0, rdm_ice(:,:)-zlicegr(:,:) ) 
     452         CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Lateral sea ice growth     [kg/m2/s] 
    458453      ENDIF 
    459454      ! 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90

    r3294 r3625  
    77 
    88   !!---------------------------------------------------------------------- 
    9    !!   lim_lat_acr_2   : lateral accretion of ice 
    10    !!---------------------------------------------------------------------- 
    11    USE par_oce          ! ocean parameters 
     9   !!   lim_lat_acr_2 : lateral accretion of ice 
     10   !!---------------------------------------------------------------------- 
     11   USE par_oce        ! ocean parameters 
    1212   USE phycst 
    1313   USE thd_ice_2 
    1414   USE ice_2 
    1515   USE limistate_2  
    16    USE lib_mpp          ! MPP library 
    17    USE wrk_nemo         ! work arrays 
     16   USE lib_mpp        ! MPP library 
     17   USE wrk_nemo       ! work arrays 
     18   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    1819 
    1920   IMPLICIT NONE 
    … …  
    145146         frld_1d   (ji) = MAX( zfrlnew , zfrlmin(ji) ) 
    146147         !--computation of the remaining part of ice thickness which has been already used 
    147          zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) ) &  
    148                       -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic )  
     148         zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) )   &  
     149            &         -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic )  
    149150      END DO 
    150151  
    … …  
    196197            &          ) / zah 
    197198          
    198          tbif_1d(ji,3) =     (  iiceform * ( zhnews2 - zdh3 )                                          * zta1  & 
     199         tbif_1d(ji,3) =     ( iiceform * ( zhnews2 - zdh3 )                                           * zta1  & 
    199200            &              + ( iiceform * zdh3 + ( 1 - iiceform ) * zdh1 )                             * zta2  & 
    200201            &              + ( iiceform * ( zhnews2 - zdh5 ) + ( 1 - iiceform ) * ( zhnews2 - zdh1 ) ) * zta3  &  
    … …  
    217218      DO ji = kideb , kiut 
    218219         dvlbq_1d  (ji) = ( 1. - frld_1d(ji) ) * h_ice_1d(ji) - ( 1. - zfrl_old(ji) ) * zhice_old(ji) 
    219          rdmicif_1d(ji) = rdmicif_1d(ji) + rhoic * dvlbq_1d(ji) 
     220         rdm_ice_1d(ji) = rdm_ice_1d(ji) + rhoic * dvlbq_1d(ji) 
     221         rdq_ice_1d(ji) = rdq_ice_1d(ji) + rcpic * dvlbq_1d(ji) * ( tfu_1d(ji) - rt0 )      ! heat content relative to rt0 
    220222      END DO 
    221223       

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90

    r3294 r3625  
    1818   USE ice_2 
    1919   USE limistate_2 
     20   USE cpl_oasis3, ONLY : lk_cpl 
    2021   USE in_out_manager 
    2122   USE lib_mpp          ! MPP library 
    2223   USE wrk_nemo         ! work arrays 
    23    USE cpl_oasis3, ONLY : lk_cpl 
    24        
     24   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     25     
    2526   IMPLICIT NONE 
    2627   PRIVATE 
    … …  
    8687      REAL(wp), POINTER, DIMENSION(:) ::   zrcpdt         ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer) 
    8788      REAL(wp), POINTER, DIMENSION(:) ::   zts_old        ! previous surface temperature 
    88       REAL(wp), POINTER, DIMENSION(:) ::   zidsn , z1midsn , zidsnic ! tempory variables 
     89      REAL(wp), POINTER, DIMENSION(:) ::   zidsn , z1midsn , zidsnic ! temporary variables 
    8990      REAL(wp), POINTER, DIMENSION(:) ::   zfnet          ! net heat flux at the top surface( incl. conductive heat flux) 
    9091      REAL(wp), POINTER, DIMENSION(:) ::   zsprecip       ! snow accumulation 
    … …  
    9899      REAL(wp), POINTER, DIMENSION(:) ::   zep            ! internal temperature of the 2nd layer of the snow/ice system 
    99100      REAL(wp), DIMENSION(3) :: &  
    100           zplediag  &    ! principle diagonal, subdiag. and supdiag. of the  
     101            zplediag  &    ! principle diagonal, subdiag. and supdiag. of the  
    101102          , zsubdiag  &    ! tri-diagonal matrix coming from the computation 
    102103          , zsupdiag  &    ! of the temperatures inside the snow-ice system 
    103104          , zsmbr          ! second member 
    104        REAL(wp) :: &  
    105           zhsu     &     ! thickness of surface layer 
    106           , zhe      &     ! effective thickness for compu. of equ. thermal conductivity 
    107           , zheshth  &     ! = zhe / thth 
    108           , zghe     &     ! correction factor of the thermal conductivity 
    109           , zumsb    &     ! parameter for numerical method to solve heat-diffusion eq. 
    110           , zkhsn    &     ! conductivity at the snow layer 
    111           , zkhic    &     ! conductivity at the ice layers 
    112           , zkint    &     ! equivalent conductivity at the snow-ice interface 
    113           , zkhsnint &     ! = zkint*dt / (hsn*rhosn*cpsn)   
    114           , zkhicint &     ! = 2*zkint*dt / (hic*rhoic*cpic) 
    115           , zpiv1 , zpiv2  &       ! tempory scalars used to solve the tri-diagonal system 
    116           , zb2 , zd2 , zb3 , zd3 & 
     105       REAL(wp) ::    & 
     106            zhsu      &    ! thickness of surface layer 
     107          , zhe       &    ! effective thickness for compu. of equ. thermal conductivity 
     108          , zheshth   &    ! = zhe / thth 
     109          , zghe      &    ! correction factor of the thermal conductivity 
     110          , zumsb     &    ! parameter for numerical method to solve heat-diffusion eq. 
     111          , zkhsn     &    ! conductivity at the snow layer 
     112          , zkhic     &    ! conductivity at the ice layers 
     113          , zkint     &    ! equivalent conductivity at the snow-ice interface 
     114          , zkhsnint  &    ! = zkint*dt / (hsn*rhosn*cpsn)   
     115          , zkhicint  &    ! = 2*zkint*dt / (hic*rhoic*cpic) 
     116          , zpiv1, zpiv2 & ! temporary scalars used to solve the tri-diagonal system 
     117          , zb2, zd2  &    ! temporary scalars used to solve the tri-diagonal system 
     118          , zb3, zd3  &    ! temporary scalars used to solve the tri-diagonal system 
    117119          , ztint          ! equivalent temperature at the snow-ice interface 
    118        REAL(wp) :: &  
    119           zexp      &     ! exponential function of the ice thickness 
    120           , zfsab     &     ! part of solar radiation stored in brine pockets 
    121           , zfts      &     ! value of energy balance function when the temp. equal surf. temp. 
    122           , zdfts     &     ! value of derivative of ztfs when the temp. equal surf. temp. 
    123           , zdts      &     ! surface temperature increment 
    124           , zqsnw_mlt &     ! energy needed to melt snow 
    125           , zdhsmlt   &     ! change in snow thickness due to melt 
    126           , zhsn      &     ! snow thickness (previous+accumulation-melt) 
    127           , zqsn_mlt_rem &  ! remaining heat coming from snow melting 
    128           , zqice_top_mlt & ! energy used to melt ice at top surface 
    129           , zdhssub      &  ! change in snow thick. due to sublimation or evaporation 
    130           , zdhisub      &  ! change in ice thick. due to sublimation or evaporation     
    131           , zdhsn        &  ! snow ice thickness increment 
    132           , zdtsn        &  ! snow internal temp. increment 
    133           , zdtic        &  ! ice internal temp. increment 
     120       REAL(wp) ::    &  
     121            zexp      &    ! exponential function of the ice thickness 
     122          , zfsab     &    ! part of solar radiation stored in brine pockets 
     123          , zfts      &    ! value of energy balance function when the temp. equal surf. temp. 
     124          , zdfts     &    ! value of derivative of ztfs when the temp. equal surf. temp. 
     125          , zdts      &    ! surface temperature increment 
     126          , zqsnw_mlt &    ! energy needed to melt snow 
     127          , zdhsmlt   &    ! change in snow thickness due to melt 
     128          , zhsn      &    ! snow thickness (previous+accumulation-melt) 
     129          , zqsn_mlt_rem & ! remaining heat coming from snow melting 
     130          , zqice_top_mlt &! energy used to melt ice at top surface 
     131          , zdhssub   &    ! change in snow thick. due to sublimation or evaporation 
     132          , zdhisub   &    ! change in ice thick. due to sublimation or evaporation     
     133          , zdhsn     &    ! snow ice thickness increment 
     134          , zdtsn     &    ! snow internal temp. increment 
     135          , zdtic     &    ! ice internal temp. increment 
    134136          , zqnes          ! conductive energy due to ice melting in the first ice layer 
    135        REAL(wp) :: &  
    136           ztbot     &      ! temperature at the bottom surface 
    137           , zfcbot    &      ! conductive heat flux at bottom surface 
    138           , zqice_bot &      ! energy used for bottom melting/growing 
    139           , zqice_bot_mlt &  ! energy used for bottom melting 
    140           , zqstbif_bot  &  ! part of energy stored in brine pockets used for bottom melting 
    141           , zqstbif_old  &  ! tempory var. for zqstbif_bot 
    142           , zdhicmlt      &  ! change in ice thickness due to bottom melting 
    143           , zdhicm        &  ! change in ice thickness var.  
    144           , zdhsnm        &  ! change in snow thickness var.  
    145           , zhsnfi        &  ! snow thickness var.  
    146           , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables 
    147           , ztb2, ztb3 
    148        REAL(wp) :: &  
    149           zdrmh         &   ! change in snow/ice thick. after snow-ice formation 
    150           , zhicnew       &   ! new ice thickness 
    151           , zhsnnew       &   ! new snow thickness 
    152           , zquot , ztneq &   ! tempory temp. variables 
    153           , zqice, zqicetot & ! total heat inside the snow/ice system 
    154           , zdfrl         &   ! change in ice concentration 
    155           , zdvsnvol      &   ! change in snow volume 
    156           , zdrfrl1, zdrfrl2 &  ! tempory scalars 
    157           , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind 
     137       REAL(wp) ::    &  
     138            ztbot     &    ! temperature at the bottom surface 
     139          , zfcbot    &    ! conductive heat flux at bottom surface 
     140          , zqice_bot &    ! energy used for bottom melting/growing 
     141          , zqice_bot_mlt &! energy used for bottom melting 
     142          , zqstbif_bot  & ! part of energy stored in brine pockets used for bottom melting 
     143          , zqstbif_old  & ! temporary var. for zqstbif_bot 
     144          , zdhicmlt  &    ! change in ice thickness due to bottom melting 
     145          , zdhicm    &    ! change in ice thickness var.  
     146          , zdhsnm    &    ! change in snow thickness var.  
     147          , zhsnfi    &    ! snow thickness var.  
     148          , zc1, zpc1 &    ! temporary variables 
     149          , zc2, zpc2 &    ! temporary variables 
     150          , zp1, zp2  &    ! temporary variables 
     151          , ztb2, ztb3     ! temporary variables 
     152       REAL(wp) ::    &  
     153            zdrmh     &    ! change in snow/ice thick. after snow-ice formation 
     154          , zhicnew   &    ! new ice thickness 
     155          , zhsnnew   &    ! new snow thickness 
     156          , zquot     & 
     157          , ztneq     &    ! temporary temp. variables 
     158          , zqice     & 
     159          , zqicetot  &    ! total heat inside the snow/ice system 
     160          , zdfrl     &    ! change in ice concentration 
     161          , zdvsnvol  &    ! change in snow volume 
     162          , zdrfrl1, zdrfrl2, zihsn, zidhb, zihic &  ! temporary scalars 
     163          , zihe, zihq, ziexp, ziqf, zihnf        &  ! temporary scalars 
     164          , zibmlt, ziqr, zihgnew, zind, ztmp        ! temporary scalars 
    158165       !!---------------------------------------------------------------------- 
    159166       CALL wrk_alloc( jpij, ztsmlt, ztbif  , zksn    , zkic    , zksndh , zfcsu  , zfcsudt , zi0      , z1mi0   , zqmax    ) 
    … …  
    169176        
    170177       DO ji = kideb , kiut 
     178          ! do nothing if the snow (ice) thickness falls below its minimum thickness 
    171179          zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) ) 
    172180          zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) ) 
    173           !--computation of energy due to surface melting 
    174           zqcmlts(ji) = ( MAX ( zzero ,  & 
    175              &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn ) 
    176           !--computation of energy due to bottom melting 
    177           zqcmltb(ji) = ( MAX( zzero , & 
    178              &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 
    179              &           + MAX( zzero , & 
    180              &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 
    181              &           ) * ( 1.0 - zihic  ) 
    182           !--limitation of  snow/ice system internal temperature 
     181          !--energy required to bring snow to its melting point (rt0_snow) 
     182          zqcmlts(ji) = ( MAX ( zzero , rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn ) 
     183          !--energy required to bring ice to its melting point (rt0_ice) 
     184          zqcmltb(ji) = ( MAX( zzero , rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )  & 
     185             &          + MAX( zzero , rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )  & 
     186             &          ) * ( 1.0 - zihic  ) 
     187          !--limitation of snow/ice system internal temperature 
    183188          tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) ) 
    184189          tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) ) 
    … …  
    480485          dvsbq_1d(ji) =  ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnw_old(ji) - zsprecip(ji) ) 
    481486          dvsbq_1d(ji) =  MIN( zzero , dvsbq_1d(ji) ) 
    482           rdmsnif_1d(ji) =  rhosn * dvsbq_1d(ji) 
     487          ztmp = rhosn * dvsbq_1d(ji) 
     488          rdm_snw_1d(ji) =  ztmp 
     489          !--heat content of the water provided to the ocean (referenced to rt0) 
     490          rdq_snw_1d(ji) =  cpic * ztmp * ( rt0_snow - rt0 ) 
    483491          !-- If the snow is completely melted the remaining heat is used to melt ice 
    484492          zqsn_mlt_rem  = MAX( zzero , -zhsn ) * xlsn 
    … …  
    623631          !---updating new ice thickness and computing the newly formed ice mass 
    624632          zhicnew   =  zihgnew * zhicnew 
    625           rdmicif_1d(ji) =  rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic 
     633          ztmp    =  ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic 
     634          rdm_ice_1d(ji) =  rdm_ice_1d(ji) + ztmp 
     635          !---heat content of the water provided to the ocean (referenced to rt0) 
     636          !   use of rt0_ice is OK for melting ice; in the case of freezing, tfu_1d should be used.  
     637          !   This is done in 9.5 section (see below) 
     638          rdq_ice_1d(ji) =  cpic * ztmp * ( rt0_ice - rt0 ) 
    626639          !---updating new snow thickness and computing the newly formed snow mass 
    627640          zhsnfi   = zhsn + zdhsnm 
    628641          h_snow_1d(ji) = MAX( zzero , zhsnfi ) 
    629           rdmsnif_1d(ji) =  rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn 
     642          ztmp = ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn 
     643          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp 
     644          !---updating the heat content of the water provided to the ocean (referenced to rt0) 
     645          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 ) 
    630646          !--remaining energy in case of total ablation 
    631647          zqocea(ji) = - ( zihsn * xlic * zdhicm + xlsn * ( zhsnfi - h_snow_1d(ji) ) ) * ( 1.0 - frld_1d(ji) ) 
    … …  
    659675          tbif_1d(ji,3) =  zihgnew * ztb3 + ( 1.0 - zihgnew ) * tfu_1d(ji) 
    660676          h_ice_1d(ji)  =  zhicnew 
     677          ! update the ice heat content given to the ocean in freezing case  
     678          ! (part due to difference between rt0_ice and tfu_1d) 
     679          ztmp = ( 1. - zidhb ) * rhoic * dvbbq_1d(ji) 
     680          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0_ice ) 
    661681       END DO 
    662682 
    … …  
    700720          dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn 
    701721          !---  volume change of ice and snow (used for ocean-ice freshwater flux computation) 
    702           rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d (ji) ) * rhoic 
    703           rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn 
     722          ztmp = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d (ji) ) * rhoic 
     723          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp 
     724          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0 ) 
     725          !!gm BUG ??   snow ==>  only needed for nn_ice_embd == 0  (standard levitating sea-ice) 
     726          ztmp = ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn 
     727          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp 
     728          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 ) 
    704729 
    705730          !---  Actualize new snow and ice thickness. 
    … …  
    748773          !--variation of ice volume and ice mass  
    749774          dvlbq_1d(ji)   = zihic * ( zfrl_old(ji) - frld_1d(ji) ) * h_ice_1d(ji) 
    750           rdmicif_1d(ji) = rdmicif_1d(ji) + dvlbq_1d(ji) * rhoic 
     775          ztmp = dvlbq_1d(ji) * rhoic 
     776          rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp 
     777!!gm 
     778!!gm   This should be split in two parts: 
     779!!gm         1-  heat required to bring sea-ice to tfu  : this part should be added to the heat flux taken from the ocean 
     780!!gm                 cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0_ice ) 
     781!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_ice_1d 
     782!!gm                 cpic * ztmp * ( rt0_ice - rt0 ) 
     783!!gm   Currently we put all the heat in rdq_ice_1d 
     784          rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0 ) 
     785          ! 
    751786          !--variation of snow volume and snow mass  
    752           zdvsnvol    = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji) 
    753           rdmsnif_1d(ji) = rdmsnif_1d(ji) + zdvsnvol * rhosn 
     787          zdvsnvol = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji) 
     788          ztmp     = zdvsnvol * rhosn 
     789          rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp 
     790!!gm 
     791!!gm   This should be split in two parts: 
     792!!gm         1-  heat required to bring snow to tfu  : this part should be added to the heat flux taken from the ocean 
     793!!gm                 cpic * ztmp * ( tbif_1d(ji,1) - rt0_snow ) 
     794!!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_snw_1d 
     795!!gm                 cpic * ztmp * ( rt0_snow - rt0 ) 
     796!!gm   Currently we put all the heat in rdq_snw_1d 
     797          rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( tbif_1d(ji,1) - rt0 ) 
     798 
    754799          h_snow_1d(ji)  = ziqf * h_snow_1d(ji) 
    755800 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90

    r3294 r3625  
    1313   !!---------------------------------------------------------------------- 
    1414   !!---------------------------------------------------------------------- 
    15    !!   lim_wri_2      : write of the diagnostics variables in ouput file  
    16    !!   lim_wri_init_2 : initialization and namelist read 
     15   !!   lim_wri_2       : write of the diagnostics variables in ouput file  
     16   !!   lim_wri_init_2  : initialization and namelist read 
    1717   !!   lim_wri_state_2 : write for initial state or/and abandon: 
    1818   !!                     > output.init.nc (if ninist = 1 in namelist) 
    … …  
    2626   USE ice_2 
    2727 
    28    USE dianam          ! build name of file (routine) 
     28   USE dianam           ! build name of file (routine) 
    2929   USE lbclnk 
    3030   USE in_out_manager 
    31    USE lib_mpp         ! MPP library 
    32    USE wrk_nemo        ! work arrays 
     31   USE lib_mpp          ! MPP library 
     32   USE wrk_nemo         ! work arrays 
    3333   USE iom 
    3434   USE ioipsl 
     35   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3536 
    3637   IMPLICIT NONE 
    … …  
    173174            zcmo(ji,jj,13) = qns(ji,jj) 
    174175            ! See thersf for the coefficient 
    175             zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ??? 
     176            zcmo(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ??? 
    176177            zcmo(ji,jj,15) = utau_ice(ji,jj) 
    177178            zcmo(ji,jj,16) = vtau_ice(ji,jj) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90

    r3294 r3625  
    118118          zcmo(ji,jj,13) = qns(ji,jj) 
    119119          ! See thersf for the coefficient 
    120           zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     120          zcmo(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
    121121          zcmo(ji,jj,15) = utau_ice(ji,jj) 
    122122          zcmo(ji,jj,16) = vtau_ice(ji,jj) 
    … …  
    161161                rcmoy(ji,jj,13) = qns(ji,jj) 
    162162                ! See thersf for the coefficient 
    163                 rcmoy(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     163                rcmoy(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
    164164                rcmoy(ji,jj,15) = utau_ice(ji,jj) 
    165165                rcmoy(ji,jj,16) = vtau_ice(ji,jj) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90

    r2715 r3625  
    6868      qstbif_1d   ,     &  !:    "                  "      qstoif 
    6969      fbif_1d     ,     &  !:    "                  "      fbif 
    70       rdmicif_1d  ,     &  !:    "                  "      rdmicif 
    71       rdmsnif_1d  ,     &  !:    "                  "      rdmsnif 
     70      rdm_ice_1d  ,     &  !:    "                  "      rdm_ice 
     71      rdq_ice_1d  ,     &  !:    "                  "      rdq_ice 
     72      rdm_snw_1d  ,     &  !:    "                  "      rdm_snw 
     73      rdq_snw_1d  ,     &  !:    "                  "      rdq_snw 
    7274      qlbbq_1d    ,     &  !:    "                  "      qlbsbq 
    7375      dmgwi_1d    ,     &  !:    "                  "      dmgwi 
    … …  
    108110         &      qstbif_1d(jpij),  fbif_1d(jpij),  Stat=ierr(2)) 
    109111         ! 
    110       ALLOCATE( rdmicif_1d(jpij), rdmsnif_1d(jpij), qlbbq_1d(jpij),   & 
     112      ALLOCATE( rdm_ice_1d(jpij), rdq_ice_1d(jpij)                  , & 
     113         &      rdm_snw_1d(jpij), rdq_snw_1d(jpij), qlbbq_1d(jpij)  , & 
    111114         &      dmgwi_1d(jpij)  , dvsbq_1d(jpij)  , rdvomif_1d(jpij), & 
    112115         &      dvbbq_1d(jpij)  , dvlbq_1d(jpij)  , dvnbq_1d(jpij)  , & 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90

    r2777 r3625  
    99   USE par_ice        ! LIM-3 parameter 
    1010   USE in_out_manager ! I/O manager 
    11    USE lib_mpp         ! MPP library 
     11   USE lib_mpp        ! MPP library 
     12   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    1213 
    1314   IMPLICIT NONE 
    … …  
    3031 
    3132   !!---------------------------------------------------------------------- 
    32    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     33   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    3334   !! $Id$ 
    3435   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r2777 r3625  
    99#if defined key_lim3 
    1010   !!---------------------------------------------------------------------- 
    11    !!   'key_lim3' :                                   LIM3 sea-ice model 
    12    !!---------------------------------------------------------------------- 
    13    USE par_ice          ! LIM sea-ice parameters 
    14    USE in_out_manager   ! I/O manager 
    15    USE lib_mpp         ! MPP library 
     11   !!   'key_lim3'                                      LIM-3 sea-ice model 
     12   !!---------------------------------------------------------------------- 
     13   USE par_ice        ! LIM sea-ice parameters 
     14   USE in_out_manager ! I/O manager 
     15   USE lib_mpp        ! MPP library 
    1616 
    1717   IMPLICIT NONE 
    … …  
    158158   !! * Share Module variables 
    159159   !!-------------------------------------------------------------------------- 
    160    INTEGER , PUBLIC ::   nstart    !: iteration number of the begining of the run  
    161    INTEGER , PUBLIC ::   nlast     !: iteration number of the end of the run  
    162    INTEGER , PUBLIC ::   nitrun    !: number of iteration 
    163    INTEGER , PUBLIC ::   numit     !: iteration number 
    164    REAL(wp), PUBLIC ::   rdt_ice   !: ice time step 
     160   INTEGER , PUBLIC ::   nstart      !: iteration number of the begining of the run  
     161   INTEGER , PUBLIC ::   nlast       !: iteration number of the end of the run  
     162   INTEGER , PUBLIC ::   nitrun      !: number of iteration 
     163   INTEGER , PUBLIC ::   numit       !: iteration number 
     164   REAL(wp), PUBLIC ::   rdt_ice     !: ice time step 
     165   REAL(wp), PUBLIC ::   r1_rdtice   !: = 1. / rdt_ice 
    165166 
    166167   !                                          !!** ice-dynamic namelist (namicedyn) ** 
    … …  
    201202   !                                              !!** ice-salinity namelist (namicesal) ** 
    202203   INTEGER , PUBLIC ::   num_sal     = 1           !: salinity configuration used in the model 
    203    !                                               ! 1 - s constant in space and time 
     204   !                                               ! 1 - constant salinity in both space and time 
    204205   !                                               ! 2 - prognostic salinity (s(z,t)) 
    205206   !                                               ! 3 - salinity profile, constant in time 
    206    !                                               ! 4 - salinity variations affect only ice thermodynamics 
    207207   INTEGER , PUBLIC ::   sal_prof    = 1           !: salinity profile or not  
    208208   INTEGER , PUBLIC ::   thcon_i_swi = 1           !: thermal conductivity: =1 Untersteiner (1964) ; =2 Pringle et al (2007) 
    … …  
    264264   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   phicif      !: Old ice thickness 
    265265   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fbif        !: Heat flux at the ice base 
    266    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmsnif     !: Variation of snow mass 
    267    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmicif     !: Variation of ice mass 
     266   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_snw     !: Variation of snow mass over 1 time step     [Kg/m2] 
     267   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_snw     !: Heat content associated with rdm_snw        [J/m2] 
     268   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_ice     !: Variation of ice mass over 1 time step      [Kg/m2] 
     269   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_ice     !: Heat content associated with rdm_ice        [J/m2] 
    268270   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qldif       !: heat balance of the lead (or of the open ocean) 
    269271   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qcmif       !: Energy needed to bring the ocean to freezing  
    … …  
    276278   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qfvbq       !: store energy in case of total lateral ablation (?) 
    277279   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dmgwi       !: Variation of the mass of snow ice 
    278    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_res   !: Residual salt flux due to correction of ice thickness 
    279    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsbri       !: Salt flux due to brine rejection 
    280    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_rpo   !: Salt flux associated with porous ridged ice formation 
    281    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_rpo   !: Heat flux associated with porous ridged ice formation 
     280   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_thd     !: salt flux due to ice growth/melt                      [PSU/m2/s] 
     281   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bri     !: salt flux due to brine rejection                      [PSU/m2/s] 
     282   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_mec     !: salt flux due to porous ridged ice formation          [PSU/m2/s] 
     283   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: residual salt flux due to correction of ice thickness [PSU/m2/s] 
    282284   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhbri       !: heat flux due to brine rejection 
    283    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmec       !: Mass flux due to snow loss during compression 
    284    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fseqv       !: Equivalent salt flux due to ice growth/melt 
    285    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhmec       !: Heat flux due to snow loss during compression 
    286    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_res   !: Residual heat flux due to correction of ice thickness 
     285   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_mec   !: heat flux associated with porous ridged ice formation [???] 
     286   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_res   !: residual heat flux due to correction of ice thickness 
     287   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmec       !: mass flux due to snow loss during compression         [Kg/m2/s] 
     288   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhmec       !: heat flux due to snow loss during compression 
    287289 
    288290   ! temporary arrays for dummy version of the code 
    … …  
    415417 
    416418   !!---------------------------------------------------------------------- 
    417    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
     419   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2010) 
    418420   !! $Id$ 
    419421   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    444446 
    445447      ii = ii + 1 
    446       ALLOCATE( firic    (jpi,jpj) , fcsic  (jpi,jpj) , fleic    (jpi,jpj) , qlatic   (jpi,jpj) ,     & 
    447          &      rdvosif  (jpi,jpj) , rdvobif(jpi,jpj) , fdvolif  (jpi,jpj) , rdvonif  (jpi,jpj) ,     & 
    448          &      sist     (jpi,jpj) , icethi (jpi,jpj) , t_bo     (jpi,jpj) , hicifp   (jpi,jpj) ,     & 
    449          &      frld     (jpi,jpj) , pfrld  (jpi,jpj) , phicif   (jpi,jpj) , fbif     (jpi,jpj) ,     & 
    450          &      rdmsnif  (jpi,jpj) , rdmicif(jpi,jpj) , qldif    (jpi,jpj) , qcmif    (jpi,jpj) ,     & 
    451          &      fdtcn    (jpi,jpj) , qdtcn  (jpi,jpj) , fstric   (jpi,jpj) , fscmbq   (jpi,jpj) ,     & 
    452          &      ffltbif  (jpi,jpj) , fsbbq  (jpi,jpj) , qfvbq    (jpi,jpj) , dmgwi    (jpi,jpj) ,     & 
    453          &      fsalt_res(jpi,jpj) , fsbri  (jpi,jpj) , fsalt_rpo(jpi,jpj) , fheat_rpo(jpi,jpj) ,     & 
    454          &      fhbri    (jpi,jpj) , fmmec  (jpi,jpj) , fseqv    (jpi,jpj) , fhmec    (jpi,jpj) ,     & 
    455          &      fheat_res(jpi,jpj)                                                              , STAT=ierr(ii) ) 
     448      ALLOCATE( firic    (jpi,jpj) , fcsic  (jpi,jpj) , fleic  (jpi,jpj) , qlatic   (jpi,jpj) ,     & 
     449         &      rdvosif  (jpi,jpj) , rdvobif(jpi,jpj) , fdvolif(jpi,jpj) , rdvonif  (jpi,jpj) ,     & 
     450         &      sist     (jpi,jpj) , icethi (jpi,jpj) , t_bo   (jpi,jpj) , hicifp   (jpi,jpj) ,     & 
     451         &      frld     (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) , fbif     (jpi,jpj) ,     & 
     452         &      rdm_snw  (jpi,jpj) , rdq_snw(jpi,jpj) , rdm_ice(jpi,jpj) , rdq_ice  (jpi,jpj) ,     & 
     453         &                                              qldif  (jpi,jpj) , qcmif    (jpi,jpj) ,     & 
     454         &      fdtcn    (jpi,jpj) , qdtcn  (jpi,jpj) , fstric (jpi,jpj) , fscmbq   (jpi,jpj) ,     & 
     455         &      ffltbif  (jpi,jpj) , fsbbq  (jpi,jpj) , qfvbq  (jpi,jpj) , dmgwi    (jpi,jpj) ,     & 
     456         &      sfx_res  (jpi,jpj) , sfx_bri(jpi,jpj) , sfx_mec(jpi,jpj) , fheat_mec(jpi,jpj) ,     & 
     457         &      fhbri    (jpi,jpj) , fmmec  (jpi,jpj) , sfx_thd(jpi,jpj) , fhmec    (jpi,jpj) ,     & 
     458         &      fheat_res(jpi,jpj)                                                            , STAT=ierr(ii) ) 
    456459 
    457460      ii = ii + 1 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90

    r3294 r3625  
    1010#if defined key_lim3 
    1111   !!---------------------------------------------------------------------- 
    12    !!   'key_lim3' :                                   LIM sea-ice model 
    13    !!---------------------------------------------------------------------- 
    14    !!   ice_init       : sea-ice model initialization 
    15    !!---------------------------------------------------------------------- 
    16    USE phycst           ! physical constants 
    17    USE dom_oce          ! ocean domain 
    18    USE sbc_oce          ! Surface boundary condition: ocean fields 
    19    USE sbc_ice          ! Surface boundary condition: ice   fields 
    20    USE ice              ! LIM variables 
    21    USE par_ice          ! LIM parameters 
    22    USE dom_ice          ! LIM domain 
    23    USE thd_ice          ! LIM thermodynamical variables 
    24    USE limitd_me        ! LIM ice thickness distribution 
    25    USE limmsh           ! LIM mesh 
    26    USE limistate        ! LIM initial state 
    27    USE limrst           ! LIM restart 
    28    USE limthd           ! LIM ice thermodynamics 
    29    USE limthd_sal       ! LIM ice thermodynamics: salinity 
    30    USE limvar           ! LIM variables 
    31    USE limsbc           ! LIM surface boundary condition 
    32    USE in_out_manager   ! I/O manager 
    33    USE lib_mpp          ! MPP library 
     12   !!   'key_lim3'                                     LIM sea-ice model 
     13   !!---------------------------------------------------------------------- 
     14   !!   ice_init      : sea-ice model initialization 
     15   !!---------------------------------------------------------------------- 
     16   USE phycst         ! physical constants 
     17   USE dom_oce        ! ocean domain 
     18   USE sbc_oce        ! Surface boundary condition: ocean fields 
     19   USE sbc_ice        ! Surface boundary condition: ice   fields 
     20   USE ice            ! LIM variables 
     21   USE par_ice        ! LIM parameters 
     22   USE dom_ice        ! LIM domain 
     23   USE thd_ice        ! LIM thermodynamical variables 
     24   USE limitd_me      ! LIM ice thickness distribution 
     25   USE limmsh         ! LIM mesh 
     26   USE limistate      ! LIM initial state 
     27   USE limrst         ! LIM restart 
     28   USE limthd         ! LIM ice thermodynamics 
     29   USE limthd_sal     ! LIM ice thermodynamics: salinity 
     30   USE limvar         ! LIM variables 
     31   USE limsbc         ! LIM surface boundary condition 
     32   USE in_out_manager ! I/O manager 
     33   USE lib_mpp        ! MPP library 
     34   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3435 
    3536   IMPLICIT NONE 
    … …  
    3940 
    4041   !!---------------------------------------------------------------------- 
    41    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     42   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    4243   !! $Id$ 
    4344   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    7980      CALL lim_thd_sal_init            ! set ice salinity parameters 
    8081      ! 
    81       rdt_ice = nn_fsbc * rdttra(1)    ! sea-ice timestep 
     82      rdt_ice   = nn_fsbc * rdttra(1)  ! sea-ice timestep 
     83      r1_rdtice = 1._wp / rdt_ice      ! sea-ice timestep inverse 
    8284      ! 
    8385      CALL lim_msh                     ! ice mesh initialization 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90

    r3294 r3625  
    1515   !!   lim_adv_y  : advection of sea ice on y axis 
    1616   !!---------------------------------------------------------------------- 
    17    USE dom_oce          ! ocean domain 
    18    USE dom_ice          ! LIM-3 domain 
    19    USE ice              ! LIM-3 variables 
    20    USE lbclnk           ! lateral boundary condition - MPP exchanges 
    21    USE in_out_manager   ! I/O manager 
    22    USE prtctl           ! Print control 
    23    USE lib_mpp          ! MPP library 
    24    USE wrk_nemo         ! work arrays 
     17   USE dom_oce        ! ocean domain 
     18   USE ice            ! LIM-3 variables 
     19   USE dom_ice        ! LIM-3 domain 
     20   USE lbclnk         ! lateral boundary condition - MPP exchanges 
     21   USE in_out_manager ! I/O manager 
     22   USE prtctl         ! Print control 
     23   USE lib_mpp        ! MPP library 
     24   USE wrk_nemo       ! work arrays 
     25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2526 
    2627   IMPLICIT NONE 
    … …  
    3738#  include "vectopt_loop_substitute.h90" 
    3839   !!---------------------------------------------------------------------- 
    39    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     40   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    4041   !! $Id$ 
    4142   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    8889            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      & 
    8990               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  ) 
    90             zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
     91            zin0    = ( 1.0 - MAX( rzero, SIGN( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
    9192 
    9293            ps0 (ji,jj) = zslpmax   
    … …  
    273274            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    274275               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  ) 
    275             zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
     276            zin0    = ( 1.0 - MAX( rzero, SIGN( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
    276277            ! 
    277278            ps0 (ji,jj) = zslpmax   

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r2777 r3625  
    1010#if defined key_lim3 
    1111   !!---------------------------------------------------------------------- 
    12    !!   'key_lim3' :                                   LIM3 sea-ice model 
     12   !!   'key_lim3'                                      LIM-3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    !!    lim_cons   :   checks whether energy, mass and salt are conserved  
     14   !!    lim_cons     :   checks whether energy, mass and salt are conserved  
    1515   !!---------------------------------------------------------------------- 
    16    USE par_ice          ! LIM-3 parameter 
    17    USE ice              ! LIM-3 variables 
    18    USE dom_ice          ! LIM-3 domain 
    19    USE dom_oce          ! ocean domain 
    20    USE in_out_manager   ! I/O manager 
    21    USE lib_mpp          ! MPP library 
     16   USE par_ice        ! LIM-3 parameter 
     17   USE ice            ! LIM-3 variables 
     18   USE dom_ice        ! LIM-3 domain 
     19   USE dom_oce        ! ocean domain 
     20   USE in_out_manager ! I/O manager 
     21   USE lib_mpp        ! MPP library 
     22   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2223 
    2324   IMPLICIT NONE 
    … …  
    2930 
    3031   !!---------------------------------------------------------------------- 
    31    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     32   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    3233   !! $Id$ 
    3334   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90

    r2715 r3625  
    1111   !!   'key_lim3'                                       LIM3 sea-ice model 
    1212   !!---------------------------------------------------------------------- 
    13    !!   lim_dia        : computation and output of the time evolution of keys variables 
    14    !!   lim_dia_init   : initialization and namelist read 
    15    !!---------------------------------------------------------------------- 
    16    USE ice             ! LIM-3: sea-ice variable 
    17    USE par_ice         ! LIM-3: ice parameters 
    18    USE dom_ice         ! LIM-3: sea-ice domain 
    19    USE dom_oce         ! ocean domain 
    20    USE sbc_oce         ! surface boundary condition: ocean fields 
    21    USE daymod          ! model calendar 
    22    USE phycst          ! physical constant 
    23    USE in_out_manager  ! I/O manager 
    24    USE lib_mpp         ! MPP library 
    25     
     13   !!   lim_dia       : computation and output of the time evolution of keys variables 
     14   !!   lim_dia_init  : initialization and namelist read 
     15   !!---------------------------------------------------------------------- 
     16   USE ice            ! LIM-3: sea-ice variable 
     17   USE par_ice        ! LIM-3: ice parameters 
     18   USE dom_ice        ! LIM-3: sea-ice domain 
     19   USE dom_oce        ! ocean domain 
     20   USE sbc_oce        ! surface boundary condition: ocean fields 
     21   USE daymod         ! model calendar 
     22   USE phycst         ! physical constant 
     23   USE in_out_manager ! I/O manager 
     24   USE lib_mpp        ! MPP library 
     25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     26 
    2627   IMPLICIT NONE 
    2728   PRIVATE 
    … …  
    7071      !!              the temporal evolution of some key variables 
    7172      !!------------------------------------------------------------------- 
    72       INTEGER  ::   jv, ji, jj, jl   ! dummy loop indices 
    73       REAL(wp) ::   zshift_date      ! date from the minimum ice extent 
    74       REAL(wp) ::   zday, zday_min   ! current day, day of minimum extent 
    75       REAL(wp) ::   zafy, zamy       ! temporary area of fy and my ice 
     73      INTEGER  ::   jv, ji, jj, jl       ! dummy loop indices 
     74      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integer 
     75      REAL(wp) ::   zshift_date          ! date from the minimum ice extent 
     76      REAL(wp) ::   zday, zday_min       ! current day, day of minimum extent 
     77      REAL(wp) ::   zafy, zamy           ! temporary area of fy and my ice 
    7678      REAL(wp) ::   zindb 
    77       REAL(wp), DIMENSION(jpinfmx) ::   vinfor           ! temporary working space  
     79      REAL(wp), DIMENSION(jpinfmx) ::   vinfor   ! 1D workspace  
    7880      !!------------------------------------------------------------------- 
    7981 
    … …  
    105107            IF( tms(ji,jj) == 1 ) THEN 
    106108               vinfor(3)  = vinfor(3)  + at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !ice area 
    107                IF (at_i(ji,jj).GT.0.15) vinfor(5) = vinfor(5) + aire(ji,jj) * 1.e-12_wp !ice extent 
     109               IF ( at_i(ji,jj) > 0.15 )  vinfor(5) = vinfor(5) + aire(ji,jj) * 1.e-12_wp !ice extent 
    108110               vinfor(7)  = vinfor(7)  + vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !ice volume 
    109111               vinfor(9)  = vinfor(9)  + vt_s(ji,jj)*aire(ji,jj) * 1.e-12_wp !snow volume 
    … …  
    111113               vinfor(29) = vinfor(29) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !mean salinity 
    112114               ! the computation of this diagnostic is not reliable 
    113                vinfor(31) = vinfor(31) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + &  
    114                   v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12  
    115                vinfor(53) = vinfor(53) + emps(ji,jj)*aire(ji,jj) * 1.e-12_wp !salt flux 
    116                vinfor(55) = vinfor(55) + fsbri(ji,jj)*aire(ji,jj) * 1.e-12_wp !brine drainage flux 
    117                vinfor(57) = vinfor(57) + fseqv(ji,jj)*aire(ji,jj) * 1.e-12_wp !equivalent salt flux 
     115               vinfor(31) = vinfor(31) + vt_i(ji,jj) * (  u_ice(ji,jj)*u_ice(ji,jj)  &  
     116                  &                                     + v_ice(ji,jj)*v_ice(ji,jj) ) * aire(ji,jj) * 1.e-12  
     117               vinfor(53) = vinfor(53) + sfx (ji,jj)*aire(ji,jj) * 1.e-12_wp !salt flux 
     118               vinfor(55) = vinfor(55) + sfx_bri(ji,jj)*aire(ji,jj) * 1.e-12_wp !brine drainage flux 
     119               vinfor(57) = vinfor(57) + sfx_thd(ji,jj)*aire(ji,jj) * 1.e-12_wp !equivalent salt flux 
    118120               vinfor(59) = vinfor(59) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SST 
    119121               vinfor(61) = vinfor(61) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SSS 
    … …  
    180182               vinfor(43) = vinfor(43) + diag_dyn_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp  
    181183               vinfor(45) = vinfor(45) + dv_dt_thd(ji,jj,5)*aire(ji,jj) * 1.e-12_wp 
    182                vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp / rdt_ice ! volume acc in OW 
     184               vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp * r1_rdtice  ! volume acc in OW 
    183185            ENDIF 
    184186         END DO 
    … …  
    231233      vinfor(51) = zindb*vinfor(51) / MAX(vinfor(27),epsi06) 
    232234 
    233       !! Fram Strait Export 
    234       !! 83 = area export 
    235       !! 84 = volume export 
    236       !! Fram strait in ORCA2 = 5 points 
    237       !! export = -v_ice*e1t*ddtb*at_i or -v_ice*e1t*ddtb*at_i*h_i 
    238       jj = 136 ! C grid 
    239       vinfor(83) = 0.0 
    240       vinfor(84) = 0.0 
    241       DO ji = 134, 138 
    242          vinfor(83) = vinfor(83) - v_ice(ji,jj) * &  
    243             e1t(ji,jj)*at_i(ji,jj)*rdt_ice * 1.e-12_wp 
    244          vinfor(84) = vinfor(84) - v_ice(ji,jj) * &  
    245             e1t(ji,jj)*vt_i(ji,jj)*rdt_ice * 1.e-12_wp 
    246       END DO 
     235      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : Fram Strait Export 
     236         SELECT CASE ( jp_cfg ) 
     237         CASE ( 2 )                          ! ORCA_R2 
     238            ij0 = 136   ;   ij1 = 136              ! Fram strait : 83 = area export 
     239            ii0 = 134   ;   ii1 = 138              !               84 = volume export 
     240            DO jj = mj0(ij0),mj1(ij1) 
     241               DO ji = mi0(ii0),mi1(ii1) 
     242                  vinfor(83) = vinfor(83) - v_ice(ji,jj) * e1t(ji,jj)*at_i(ji,jj)*rdt_ice * 1.e-12_wp 
     243                  vinfor(84) = vinfor(84) - v_ice(ji,jj) * e1t(ji,jj)*vt_i(ji,jj)*rdt_ice * 1.e-12_wp 
     244               END DO 
     245            END DO 
     246         END SELECT 
     247!!gm   just above, this is NOT the correct way of evaluating the transport ! 
     248!!gm        mass of snow is missing and v_ice should be the mean between jj and jj+1 
     249!!gm   Other ORCA configurations should be added 
     250      ENDIF 
    247251 
    248252      !!------------------------------------------------------------------- 
    … …  
    264268               vinfor(32) = vinfor(32) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + &  
    265269                  v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel 
    266                vinfor(54) = vinfor(54) + at_i(ji,jj)*emps(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Total salt flux 
    267                vinfor(56) = vinfor(56) + at_i(ji,jj)*fsbri(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Brine drainage salt flux 
    268                vinfor(58) = vinfor(58) + at_i(ji,jj)*fseqv(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Equivalent salt flux 
     270               vinfor(54) = vinfor(54) + sfx (ji,jj)*aire(ji,jj) * 1.e-12_wp ! Total salt flux 
     271               vinfor(56) = vinfor(56) + sfx_bri(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Brine drainage salt flux 
     272               vinfor(58) = vinfor(58) + sfx_thd(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Equivalent salt flux 
    269273               vinfor(60) = vinfor(60) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SST 
    270274               vinfor(62) = vinfor(62) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SSS 
    … …  
    331335               vinfor(44) = vinfor(44) + diag_dyn_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp  
    332336               vinfor(46) = vinfor(46) + dv_dt_thd(ji,jj,5)*aire(ji,jj) * 1.e-12_wp 
    333                vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp / rdt_ice ! volume acc in OW 
     337               vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp * r1_rdtice  ! volume acc in OW 
    334338            ENDIF 
    335339         END DO 
    … …  
    345349         END DO 
    346350      END DO 
    347       zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4))) ! 
    348       vinfor(64) = zindb * vinfor(64) / MAX(vinfor(4),epsi06) ! divide by ice extt 
     351      zindb      = 1._wp - MAX(  0._wp , SIGN( 1._wp , -vinfor(4) )  )  ! 
     352      vinfor(64) = zindb * vinfor(64) / MAX( vinfor(4) , epsi06 )  ! divide by ice extt 
    349353      !! 2.2) Diagnostics dependent on age 
    350354      !!------------------------------------ 
    … …  
    368372                  ENDIF 
    369373               END DO ! jl 
    370                IF ((at_i(ji,jj).GT.0.15).AND.(zafy.GT.zamy)) THEN 
     374               IF ( at_i(ji,jj)  >  0.15  .AND. zafy  >  zamy ) THEN 
    371375                  vinfor(22) = vinfor(22) + aire(ji,jj) * 1.e-12_wp ! Seasonal ice extent 
    372376               ENDIF 
    373                IF ((at_i(ji,jj).GT.0.15).AND.(zafy.LE.zamy)) THEN 
     377               IF ( at_i(ji,jj)  >  0.15  .AND. zafy <= zamy ) THEN 
    374378                  vinfor(24) = vinfor(24) + aire(ji,jj) * 1.e-12_wp ! Perennial ice extent 
    375379               ENDIF 
    … …  
    377381         END DO ! jj 
    378382      END DO ! ji 
    379       zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(26))) !=0 if no multiyear ice 1 if yes 
    380       vinfor(50) = zindb*vinfor(50) / MAX(vinfor(26),epsi06) 
    381       zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(28))) !=0 if no multiyear ice 1 if yes 
    382       vinfor(52) = zindb*vinfor(52) / MAX(vinfor(28),epsi06) 
     383      zindb      = 1.0 - MAX(  0.0,SIGN( 1._wp , -vinfor(26) )  )    !=0 if no multiyear ice 1 if yes 
     384      vinfor(50) = zindb * vinfor(50) / MAX( vinfor(26) , epsi06 ) 
     385      zindb      = 1.0 - MAX(  0._wp , SIGN( 1._wp , -vinfor(28) )  ) !=0 if no multiyear ice 1 if yes 
     386      vinfor(52) = zindb * vinfor(52) / MAX( vinfor(28) , epsi06 ) 
    383387 
    384388      !  Accumulation before averaging  

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r3294 r3625  
    1515   !!    lim_dyn_init : initialization and namelist read 
    1616   !!---------------------------------------------------------------------- 
    17    USE phycst           ! physical constants 
    18    USE dom_oce          ! ocean space and time domain 
    19    USE sbc_oce          ! Surface boundary condition: ocean fields 
    20    USE sbc_ice          ! Surface boundary condition: ice   fields 
    21    USE ice              ! LIM-3 variables 
    22    USE par_ice          ! LIM-3 parameters 
    23    USE dom_ice          ! LIM-3 domain 
    24    USE limrhg           ! LIM-3 rheology 
    25    USE lbclnk           ! lateral boundary conditions - MPP exchanges 
    26    USE lib_mpp          ! MPP library 
    27    USE wrk_nemo         ! work arrays 
    28    USE in_out_manager   ! I/O manager 
    29    USE prtctl           ! Print control 
     17   USE phycst         ! physical constants 
     18   USE dom_oce        ! ocean space and time domain 
     19   USE sbc_oce        ! Surface boundary condition: ocean fields 
     20   USE sbc_ice        ! Surface boundary condition: ice   fields 
     21   USE ice            ! LIM-3 variables 
     22   USE par_ice        ! LIM-3 parameters 
     23   USE dom_ice        ! LIM-3 domain 
     24   USE limrhg         ! LIM-3 rheology 
     25   USE lbclnk         ! lateral boundary conditions - MPP exchanges 
     26   USE lib_mpp        ! MPP library 
     27   USE wrk_nemo       ! work arrays 
     28   USE in_out_manager ! I/O manager 
     29   USE prtctl         ! Print control 
     30   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3031 
    3132   IMPLICIT NONE 
    … …  
    3738#  include "vectopt_loop_substitute.h90" 
    3839   !!---------------------------------------------------------------------- 
    39    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     40   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    4041   !! $Id$ 
    4142   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r3294 r3625  
    1212   !!   'key_lim3'                                      LIM3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    !!   lim_hdf  : diffusion trend on sea-ice variable 
     14   !!   lim_hdf       : diffusion trend on sea-ice variable 
    1515   !!---------------------------------------------------------------------- 
    16    USE dom_oce          ! ocean domain 
    17    USE ice              ! LIM-3: ice variables 
    18    USE lbclnk           ! lateral boundary condition - MPP exchanges 
    19    USE lib_mpp          ! MPP library 
    20    USE wrk_nemo         ! work arrays 
    21    USE prtctl           ! Print control 
    22    USE in_out_manager   ! I/O manager 
     16   USE dom_oce        ! ocean domain 
     17   USE ice            ! LIM-3: ice variables 
     18   USE lbclnk         ! lateral boundary condition - MPP exchanges 
     19   USE lib_mpp        ! MPP library 
     20   USE wrk_nemo       ! work arrays 
     21   USE prtctl         ! Print control 
     22   USE in_out_manager ! I/O manager 
     23   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2324 
    2425   IMPLICIT NONE 
    … …  
    3435#  include "vectopt_loop_substitute.h90" 
    3536   !!---------------------------------------------------------------------- 
    36    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
     37   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2010) 
    3738   !! $Id$ 
    3839   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r3610 r3625  
    2626   USE lib_mpp          ! MPP library 
    2727   USE wrk_nemo         ! work arrays 
     28   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2829 
    2930   IMPLICIT NONE 
    … …  
    4849 
    4950   !!---------------------------------------------------------------------- 
    50    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     51   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    5152   !! $Id$ 
    5253   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r3294 r3625  
    1010#if defined key_lim3 
    1111   !!---------------------------------------------------------------------- 
    12    !!   'key_lim3' :                                    LIM3 sea-ice model 
     12   !!   'key_lim3'                                      LIM-3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    USE par_oce          ! ocean parameters 
    15    USE dom_oce          ! ocean domain 
    16    USE phycst           ! physical constants (ocean directory)  
    17    USE sbc_oce          ! surface boundary condition: ocean fields 
    18    USE thd_ice          ! LIM thermodynamics 
    19    USE ice              ! LIM variables 
    20    USE par_ice          ! LIM parameters 
    21    USE dom_ice          ! LIM domain 
    22    USE limthd_lac       ! LIM 
    23    USE limvar           ! LIM 
    24    USE limcons          ! LIM 
    25    USE in_out_manager   ! I/O manager 
    26    USE lbclnk           ! lateral boundary condition - MPP exchanges 
    27    USE lib_mpp          ! MPP library 
    28    USE wrk_nemo         ! work arrays 
    29    USE prtctl           ! Print control 
     14   USE par_oce        ! ocean parameters 
     15   USE dom_oce        ! ocean domain 
     16   USE phycst         ! physical constants (ocean directory)  
     17   USE sbc_oce        ! surface boundary condition: ocean fields 
     18   USE thd_ice        ! LIM thermodynamics 
     19   USE ice            ! LIM variables 
     20   USE par_ice        ! LIM parameters 
     21   USE dom_ice        ! LIM domain 
     22   USE limthd_lac     ! LIM 
     23   USE limvar         ! LIM 
     24   USE limcons        ! LIM 
     25   USE in_out_manager ! I/O manager 
     26   USE lbclnk         ! lateral boundary condition - MPP exchanges 
     27   USE lib_mpp        ! MPP library 
     28   USE wrk_nemo       ! work arrays 
     29   USE prtctl         ! Print control 
     30   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3031 
    3132   IMPLICIT NONE 
    … …  
    3839   PUBLIC   lim_itd_me_alloc        ! called by iceini.F90 
    3940 
    40    REAL(wp)  ::   epsi11 = 1.e-11_wp   ! constant values 
    41    REAL(wp)  ::   epsi10 = 1.e-10_wp   ! constant values 
    42    REAL(wp)  ::   epsi06 = 1.e-06_wp   ! constant values 
     41   REAL(wp) ::   epsi11 = 1.e-11_wp   ! constant values 
     42   REAL(wp) ::   epsi10 = 1.e-10_wp   ! constant values 
     43   REAL(wp) ::   epsi06 = 1.e-06_wp   ! constant values 
    4344 
    4445   !----------------------------------------------------------------------- 
    … …  
    4748   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
    4849   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
    49  
     50   ! 
    5051   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/ 
    5152   !                                                           !  closing associated w/ category n 
    52  
     53   ! 
    5354   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
    5455   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
    … …  
    7071   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dvirdgdt   ! rate of ice volume ridged (m/s) 
    7172   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   opening    ! rate of opening due to divergence/shear (1/s) 
     73   ! 
    7274   !!---------------------------------------------------------------------- 
    7375   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
    … …  
    126128      INTEGER ::   ji, jj, jk, jl   ! dummy loop index 
    127129      INTEGER ::   niter, nitermax = 20   ! local integer  
    128       LOGICAL  ::   asum_error              ! flag for asum .ne. 1 
    129       INTEGER  ::   iterate_ridging         ! if true, repeat the ridging 
    130       REAL(wp) ::   w1, tmpfac, dti         ! local scalar 
     130      LOGICAL  ::   asum_error            ! flag for asum .ne. 1 
     131      INTEGER  ::   iterate_ridging       ! if true, repeat the ridging 
     132      REAL(wp) ::   w1, tmpfac            ! local scalar 
    131133      CHARACTER (len = 15) ::   fieldid 
    132134      REAL(wp), POINTER, DIMENSION(:,:) ::   closing_net     ! net rate at which area is removed    (1/s) 
    … …  
    152154      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    153155      !-----------------------------------------------------------------------------! 
    154       ! Set hi_max(ncat) to a big value to ensure that all ridged ice  
    155       ! is thinner than hi_max(ncat). 
     156      ! Set hi_max(ncat) to a big value to ensure that all ridged ice is thinner than hi_max(ncat). 
    156157 
    157158      hi_max(jpl) = 999.99 
    158159 
    159       Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0      ! proport const for PE 
    160       CALL lim_itd_me_ridgeprep ! prepare ridging 
    161  
     160      Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0                ! proport const for PE 
     161      ! 
     162      CALL lim_itd_me_ridgeprep                                    ! prepare ridging 
     163      ! 
    162164      IF( con_i)   CALL lim_column_sum( jpl, v_i, vt_i_init )      ! conservation check 
    163165 
    … …  
    166168            msnow_mlt(ji,jj) = 0._wp 
    167169            esnow_mlt(ji,jj) = 0._wp 
    168             dardg1dt (ji,jj)  = 0._wp 
    169             dardg2dt (ji,jj)  = 0._wp 
    170             dvirdgdt (ji,jj)  = 0._wp 
    171             opening  (ji,jj)  = 0._wp 
     170            dardg1dt (ji,jj) = 0._wp 
     171            dardg2dt (ji,jj) = 0._wp 
     172            dvirdgdt (ji,jj) = 0._wp 
     173            opening  (ji,jj) = 0._wp 
    172174 
    173175            !-----------------------------------------------------------------------------! 
    … …  
    201203            ! to give asum = 1.0 after ridging. 
    202204 
    203             divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) / rdt_ice  ! asum found in ridgeprep 
     205            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep 
    204206 
    205207            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) 
    … …  
    207209            ! 2.3 opning 
    208210            !------------ 
    209             ! Compute the (non-negative) opening rate that will give  
    210             ! asum = 1.0 after ridging. 
     211            ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging. 
    211212            opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj) 
    212213         END DO 
    … …  
    257258                  IF ( a_i(ji,jj,jl) > epsi11 .AND. athorn(ji,jj,jl) > 0._wp )THEN 
    258259                     w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
    259                      IF ( w1 > a_i(ji,jj,jl) ) THEN 
     260                     IF ( w1  > a_i(ji,jj,jl) ) THEN 
    260261                        tmpfac = a_i(ji,jj,jl) / w1 
    261262                        closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 
    … …  
    291292               ELSE 
    292293                  iterate_ridging    = 1 
    293                   divu_adv   (ji,jj) = (1._wp - asum(ji,jj)) / rdt_ice 
     294                  divu_adv   (ji,jj) = (1._wp - asum(ji,jj)) * r1_rdtice 
    294295                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 
    295296                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) ) 
    … …  
    308309 
    309310         IF( iterate_ridging == 1 ) THEN 
    310             IF( niter .GT. nitermax ) THEN 
     311            IF( niter  > nitermax ) THEN 
    311312               WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 
    312313               WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging 
    … …  
    323324      ! Update fresh water and heat fluxes due to snow melt. 
    324325 
    325       dti = 1._wp / rdt_ice 
    326  
    327326      asum_error = .false.  
    328327 
    … …  
    330329         DO ji = 1, jpi 
    331330 
    332             IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true. 
    333  
    334             dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti 
    335             dardg2dt(ji,jj) = dardg2dt(ji,jj) * dti 
    336             dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * dti 
    337             opening (ji,jj) = opening (ji,jj) * dti 
     331            IF( ABS( asum(ji,jj) - 1.0 ) > epsi11 )  asum_error = .true. 
     332 
     333            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice 
     334            dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice 
     335            dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice 
     336            opening (ji,jj) = opening (ji,jj) * r1_rdtice 
    338337 
    339338            !-----------------------------------------------------------------------------! 
    340339            ! 5) Heat, salt and freshwater fluxes 
    341340            !-----------------------------------------------------------------------------! 
    342             fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * dti     ! fresh water source for ocean 
    343             fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * dti     ! heat sink for ocean 
     341            fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
     342            fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean 
    344343 
    345344         END DO 
    … …  
    349348      DO jj = 1, jpj 
    350349         DO ji = 1, jpi 
    351             IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) THEN ! there is a bug 
     350            IF( ABS( asum(ji,jj) - 1._wp )  >  epsi11 ) THEN  ! there is a bug 
    352351               WRITE(numout,*) ' ' 
    353352               WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
    … …  
    391390      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
    392391      d_smv_i_trp(:,:,:)   = 0._wp 
    393       IF(  num_sal == 2  .OR.  num_sal == 4  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
     392      IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
    394393 
    395394      IF(ln_ctl) THEN     ! Control print 
    … …  
    430429 
    431430      ! update of fields will be made later in lim update 
    432       u_ice(:,:)    = old_u_ice(:,:) 
    433       v_ice(:,:)    = old_v_ice(:,:) 
    434       a_i(:,:,:)    = old_a_i(:,:,:) 
    435       v_s(:,:,:)    = old_v_s(:,:,:) 
    436       v_i(:,:,:)    = old_v_i(:,:,:) 
    437       e_s(:,:,:,:)  = old_e_s(:,:,:,:) 
    438       e_i(:,:,:,:)  = old_e_i(:,:,:,:) 
    439       oa_i(:,:,:)   = old_oa_i(:,:,:) 
    440       IF(  num_sal == 2  .OR.  num_sal == 4  )   smv_i(:,:,:) = old_smv_i(:,:,:) 
     431      u_ice(:,:)   = old_u_ice(:,:) 
     432      v_ice(:,:)   = old_v_ice(:,:) 
     433      a_i(:,:,:)   = old_a_i(:,:,:) 
     434      v_s(:,:,:)   = old_v_s(:,:,:) 
     435      v_i(:,:,:)   = old_v_i(:,:,:) 
     436      e_s(:,:,:,:) = old_e_s(:,:,:,:) 
     437      e_i(:,:,:,:) = old_e_i(:,:,:,:) 
     438      oa_i(:,:,:)  = old_oa_i(:,:,:) 
     439      IF(  num_sal == 2  )   smv_i(:,:,:) = old_smv_i(:,:,:) 
    441440 
    442441      !----------------------------------------------------! 
    … …  
    465464         DO jj = 1, jpj 
    466465            DO ji = 1, jpi 
    467                IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 
    468                   ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 
    469                   old_v_i(ji,jj,jl)     = d_v_i_trp(ji,jj,jl) 
    470                   d_v_i_trp(ji,jj,jl)   = 0._wp 
    471                   old_a_i(ji,jj,jl)     = d_a_i_trp(ji,jj,jl) 
    472                   d_a_i_trp(ji,jj,jl)   = 0._wp 
    473                   old_v_s(ji,jj,jl)     = d_v_s_trp(ji,jj,jl) 
    474                   d_v_s_trp(ji,jj,jl)   = 0._wp 
    475                   old_e_s(ji,jj,1,jl)  = d_e_s_trp(ji,jj,1,jl) 
    476                   d_e_s_trp(ji,jj,1,jl) = 0._wp 
    477                   old_oa_i(ji,jj,jl)    = d_oa_i_trp(ji,jj,jl) 
    478                   d_oa_i_trp(ji,jj,jl)  = 0._wp 
    479                   IF(  num_sal == 2  .OR.  num_sal == 4  )   old_smv_i(ji,jj,jl)  = d_smv_i_trp(ji,jj,jl) 
    480                   d_smv_i_trp(ji,jj,jl) = 0._wp 
     466               IF(  old_v_i  (ji,jj,jl) < epsi06 .AND. & 
     467                    d_v_i_trp(ji,jj,jl) > epsi06    ) THEN 
     468                  old_v_i   (ji,jj,jl)   = d_v_i_trp(ji,jj,jl) 
     469                  d_v_i_trp (ji,jj,jl)   = 0._wp 
     470                  old_a_i   (ji,jj,jl)   = d_a_i_trp(ji,jj,jl) 
     471                  d_a_i_trp (ji,jj,jl)   = 0._wp 
     472                  old_v_s   (ji,jj,jl)   = d_v_s_trp(ji,jj,jl) 
     473                  d_v_s_trp (ji,jj,jl)   = 0._wp 
     474                  old_e_s   (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 
     475                  d_e_s_trp (ji,jj,1,jl) = 0._wp 
     476                  old_oa_i  (ji,jj,jl)   = d_oa_i_trp(ji,jj,jl) 
     477                  d_oa_i_trp(ji,jj,jl)   = 0._wp 
     478                  IF(  num_sal == 2  )   old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 
     479                  d_smv_i_trp(ji,jj,jl)  = 0._wp 
    481480               ENDIF 
    482481            END DO 
    483482         END DO 
    484483      END DO 
    485  
     484      ! 
    486485      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    487486      ! 
    … …  
    612611                  ! present 
    613612                  zworka(ji,jj) = 4.0 * strength(ji,jj)              & 
    614                      + strength(ji-1,jj) * tms(ji-1,jj) &   
    615                      + strength(ji+1,jj) * tms(ji+1,jj) &   
    616                      + strength(ji,jj-1) * tms(ji,jj-1) &   
    617                      + strength(ji,jj+1) * tms(ji,jj+1)     
     613                     &          + strength(ji-1,jj) * tms(ji-1,jj) &   
     614                     &          + strength(ji+1,jj) * tms(ji+1,jj) &   
     615                     &          + strength(ji,jj-1) * tms(ji,jj-1) &   
     616                     &          + strength(ji,jj+1) * tms(ji,jj+1)     
    618617 
    619618                  zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1) 
    620619                  zworka(ji,jj) = zworka(ji,jj) / zw1 
    621620               ELSE 
    622                   zworka(ji,jj) = 0.0 
     621                  zworka(ji,jj) = 0._wp 
    623622               ENDIF 
    624623            END DO 
    … …  
    10481047         DO jj = 1, jpj 
    10491048            DO ji = 1, jpi 
    1050                IF (aicen_init(ji,jj,jl1) .GT. epsi11 .AND. athorn(ji,jj,jl1) .GT. 0.0       & 
    1051                   .AND. closing_gross(ji,jj) > 0.0) THEN 
     1049               IF( aicen_init(ji,jj,jl1)  >  epsi11  .AND.  athorn(ji,jj,jl1)  >  0._wp       & 
     1050                  .AND. closing_gross(ji,jj) > 0._wp ) THEN 
    10521051                  icells = icells + 1 
    10531052                  indxi(icells) = ji 
    … …  
    11301129            ! Salinity 
    11311130            !------------- 
    1132             smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of water frozen in voids 
     1131            smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of seawater frozen in voids 
    11331132 
    11341133            zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
    … …  
    11371136             
    11381137            !                                                             ! excess of salt is flushed into the ocean 
    1139             fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic / rdt_ice 
    1140  
     1138            sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 
     1139 
     1140            rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic / rau0   ! increase in ice volume du to seawater frozen in voids 
     1141             
    11411142            !------------------------------------             
    11421143            ! 3.6 Increment ridging diagnostics 
    … …  
    11481149            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 
    11491150            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 
    1150             diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) / rdt_ice 
    1151             opening    (ji,jj) = opening (ji,jj) + opning(ji,jj)*rdt_ice 
     1151            diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 
     1152            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 
    11521153 
    11531154            IF( con_i )   vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj) 
    … …  
    11561157            ! 3.7 Put the snow somewhere in the ocean 
    11571158            !------------------------------------------             
    1158  
    11591159            !  Place part of the snow lost by ridging into the ocean.  
    11601160            !  Note that esnow_mlt < 0; the ocean must cool to melt snow. 
    … …  
    11791179            !           ij looping 1-icells 
    11801180 
    1181             dhr(ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) 
     1181            dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) 
    11821182            dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 
    1183  
    11841183 
    11851184         END DO                 ! ij 
    … …  
    12111210 
    12121211               ! heat flux 
    1213                fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / rdt_ice 
     1212               fheat_mec(ji,jj) = fheat_mec(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) * r1_rdtice 
    12141213 
    12151214               ! Correct dimensions to avoid big values 
    … …  
    12751274               ! Transfer area, volume, and energy accordingly. 
    12761275 
    1277                IF (hrmin(ji,jj,jl1) .GE. hi_max(jl2) .OR.        & 
    1278                   hrmax(ji,jj,jl1) .LE. hi_max(jl2-1)) THEN 
    1279                   hL = 0.0 
    1280                   hR = 0.0 
     1276               IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR.        & 
     1277                   hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 
     1278                  hL = 0._wp 
     1279                  hR = 0._wp 
    12811280               ELSE 
    1282                   hL = MAX (hrmin(ji,jj,jl1), hi_max(jl2-1)) 
    1283                   hR = MIN (hrmax(ji,jj,jl1), hi_max(jl2)) 
     1281                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) ) 
     1282                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   ) 
    12841283               ENDIF 
    12851284 
    12861285               ! fraction of ridged ice area and volume going to n2 
    1287                farea = (hR-hL) / dhr(ji,jj)  
    1288                fvol(ji,jj) = (hR*hR - hL*hL) / dhr2(ji,jj) 
    1289  
    1290                a_i  (ji,jj,jl2)   = a_i  (ji,jj,jl2)  + ardg2 (ji,jj) * farea 
    1291                v_i  (ji,jj,jl2)   = v_i  (ji,jj,jl2)  + vrdg2 (ji,jj) * fvol(ji,jj) 
    1292                v_s  (ji,jj,jl2)   = v_s  (ji,jj,jl2)  + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg 
     1286               farea = ( hR - hL ) / dhr(ji,jj)  
     1287               fvol(ji,jj) = ( hR*hR - hL*hL ) / dhr2(ji,jj) 
     1288 
     1289               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ardg2 (ji,jj) * farea 
     1290               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj) 
     1291               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg 
    12931292               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg 
    1294                smv_i(ji,jj,jl2)   = smv_i(ji,jj,jl2)  + srdg2 (ji,jj) * fvol(ji,jj) 
    1295                oa_i (ji,jj,jl2)   = oa_i (ji,jj,jl2)  + oirdg2(ji,jj) * farea 
     1293               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + srdg2 (ji,jj) * fvol(ji,jj) 
     1294               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea 
    12961295 
    12971296            END DO ! ij 
    … …  
    13171316               ! Compute the fraction of rafted ice area and volume going to  
    13181317               ! thickness category jl2, transfer area, volume, and energy accordingly. 
    1319  
    1320                IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND.        & 
    1321                   hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 
    1322                   a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + arft2(ji,jj) 
    1323                   v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + virft(ji,jj) 
    1324                   v_s(ji,jj,jl2) = v_s(ji,jj,jl2) + vsrft(ji,jj)*fsnowrft 
    1325                   e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) + esrft(ji,jj)*fsnowrft 
    1326                   smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + smrft(ji,jj)     
    1327                   oa_i(ji,jj,jl2)  = oa_i(ji,jj,jl2) + oirft2(ji,jj)     
     1318               ! 
     1319               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND.        & 
     1320                   hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
     1321                  a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj) 
     1322                  v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj) 
     1323                  v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * fsnowrft 
     1324                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft 
     1325                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)     
     1326                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)     
    13281327               ENDIF ! hraft 
    1329  
     1328               ! 
    13301329            END DO ! ij 
    13311330 
    … …  
    13361335                  ji = indxi(ij) 
    13371336                  jj = indxj(ij) 
    1338                   IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND.        & 
    1339                      hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 
     1337                  IF(  hraft(ji,jj,jl1)  <=  hi_max(jl2)  .AND.        & 
     1338                       hraft(ji,jj,jl1)  >   hi_max(jl2-1)  ) THEN 
    13401339                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 
    13411340                  ENDIF 
    … …  
    15041503            DO jj = 1 , jpj 
    15051504               DO ji = 1 , jpi 
    1506 !!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice 
     1505!!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice 
    15071506!!gm                  xtmp = xtmp * unit_fac 
    15081507                  ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
    … …  
    15241523               ! fluxes are positive to the ocean 
    15251524               ! here the flux has to be negative for the ocean 
    1526 !!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice 
     1525!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice 
    15271526               !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
    15281527 
    1529 !!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice !RB   ??????? 
     1528!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ??????? 
    15301529 
    15311530               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) ) 
    … …  
    15361535 
    15371536               !           xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt 
    1538                !           fresh(i,j)      = fresh(i,j)      + xtmp 
    1539                !           fresh_hist(i,j) = fresh_hist(i,j) + xtmp 
    1540  
    1541                !           fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj)                  ) * &  
    1542                !                               rhosn * v_s(ji,jj,jl) / rdt_ice 
    1543  
    1544                !           fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * &  
    1545                !                               rhoic * v_i(ji,jj,jl) / rdt_ice 
    1546  
    1547                !           emps(i,j)      = emps(i,j)      + xtmp 
    1548                !           fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp 
     1537               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj)                  )   & 
     1538               !                                            * rhosn * v_s(ji,jj,jl) * r1_rdtice 
     1539               !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) )   &  
     1540               !                                            * rhoic * v_i(ji,jj,jl) * r1_rdtice 
     1541               !           sfx (i,j)      = sfx (i,j)      + xtmp 
    15491542 
    15501543               ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r3294 r3625  
    22   !!====================================================================== 
    33   !!                       ***  MODULE limitd_th *** 
    4    !!              Thermodynamics of ice thickness distribution 
    5    !!                   computation of changes in g(h)       
     4   !!   LIM3 ice model : ice thickness distribution: Thermodynamics 
    65   !!====================================================================== 
    76   !! History :   -   !          (W. H. Lipscomb and E.C. Hunke) CICE (c) original code 
    … …  
    2019   !!   lim_itd_shiftice : 
    2120   !!---------------------------------------------------------------------- 
    22    USE dom_ice          ! LIM-3 domain 
    23    USE par_oce          ! ocean parameters 
    24    USE dom_oce          ! ocean domain 
    25    USE phycst           ! physical constants (ocean directory)  
    26    USE thd_ice          ! LIM-3 thermodynamic variables 
    27    USE ice              ! LIM-3 variables 
    28    USE par_ice          ! LIM-3 parameters 
    29    USE limthd_lac       ! LIM-3 lateral accretion 
    30    USE limvar           ! LIM-3 variables 
    31    USE limcons          ! LIM-3 conservation 
    32    USE prtctl           ! Print control 
    33    USE in_out_manager   ! I/O manager 
    34    USE lib_mpp          ! MPP library 
    35    USE wrk_nemo         ! work arrays 
     21   USE par_oce        ! ocean parameters 
     22   USE dom_oce        ! ocean domain 
     23   USE phycst         ! physical constants (ocean directory)  
     24   USE ice            ! LIM-3 variables 
     25   USE par_ice        ! LIM-3 parameters 
     26   USE dom_ice        ! LIM-3 domain 
     27   USE thd_ice        ! LIM-3 thermodynamic variables 
     28   USE limthd_lac     ! LIM-3 lateral accretion 
     29   USE limvar         ! LIM-3 variables 
     30   USE limcons        ! LIM-3 conservation 
     31   USE prtctl         ! Print control 
     32   USE in_out_manager ! I/O manager 
     33   USE lib_mpp        ! MPP library 
     34   USE wrk_nemo       ! work arrays 
     35   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3636 
    3737   IMPLICIT NONE 
    … …  
    4949 
    5050   !!---------------------------------------------------------------------- 
    51    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
     51   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2010) 
    5252   !! $Id$ 
    5353   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    101101 
    102102      !- Trend terms 
    103       d_a_i_thd (:,:,:)  = a_i(:,:,:)   - old_a_i(:,:,:)  
    104       d_v_s_thd (:,:,:)  = v_s(:,:,:)   - old_v_s(:,:,:) 
    105       d_v_i_thd (:,:,:)  = v_i(:,:,:)   - old_v_i(:,:,:)   
     103      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:)  
     104      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:) 
     105      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)   
    106106      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)  
    107107      d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 
    108108 
    109109      d_smv_i_thd(:,:,:) = 0._wp 
    110       IF( num_sal == 2 .OR. num_sal == 4 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     110      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
    111111 
    112112      IF(ln_ctl) THEN   ! Control print 
    … …  
    143143 
    144144      !- Recover Old values 
    145       a_i(:,:,:)   = old_a_i (:,:,:) 
    146       v_s(:,:,:)   = old_v_s (:,:,:) 
    147       v_i(:,:,:)   = old_v_i (:,:,:) 
    148       e_s(:,:,:,:) = old_e_s (:,:,:,:) 
    149       e_i(:,:,:,:) = old_e_i (:,:,:,:) 
    150       ! 
    151       IF( num_sal == 2 .OR. num_sal == 4 )   smv_i(:,:,:)       = old_smv_i (:,:,:) 
     145      a_i(:,:,:)   = old_a_i(:,:,:) 
     146      v_s(:,:,:)   = old_v_s(:,:,:) 
     147      v_i(:,:,:)   = old_v_i(:,:,:) 
     148      e_s(:,:,:,:) = old_e_s(:,:,:,:) 
     149      e_i(:,:,:,:) = old_e_i(:,:,:,:) 
     150      ! 
     151      IF( num_sal == 2 )   smv_i(:,:,:) = old_smv_i(:,:,:) 
    152152      ! 
    153153   END SUBROUTINE lim_itd_th 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90

    r2715 r3625  
    1010   !!   'key_lim3'                                      LIM3 sea-ice model 
    1111   !!---------------------------------------------------------------------- 
    12    !!   lim_msh   : definition of the ice mesh 
     12   !!   lim_msh       : definition of the ice mesh 
    1313   !!---------------------------------------------------------------------- 
    1414   USE phycst         ! physical constants 
    … …  
    1818   USE lbclnk         ! lateral boundary condition - MPP exchanges 
    1919   USE lib_mpp        ! MPP library 
     20   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2021 
    2122   IMPLICIT NONE 
    … …  
    2526 
    2627   !!---------------------------------------------------------------------- 
    27    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     28   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    2829   !! $Id$ 
    2930   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r3294 r3625  
    1515   !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model 
    1616   !!---------------------------------------------------------------------- 
    17    !!   lim_rhg   : computes ice velocities 
     17   !!   lim_rhg       : computes ice velocities 
    1818   !!---------------------------------------------------------------------- 
    19    USE phycst           ! Physical constant 
    20    USE par_oce          ! Ocean parameters 
    21    USE dom_oce          ! Ocean domain 
    22    USE sbc_oce          ! Surface boundary condition: ocean fields 
    23    USE sbc_ice          ! Surface boundary condition: ice fields 
    24    USE lbclnk           ! Lateral Boundary Condition / MPP link 
    25    USE lib_mpp          ! MPP library 
    26    USE wrk_nemo         ! work arrays 
    27    USE in_out_manager   ! I/O manager 
    28    USE prtctl           ! Print control 
     19   USE phycst         ! Physical constant 
     20   USE oce     , ONLY :  snwice_mass, snwice_mass_b 
     21   USE par_oce        ! Ocean parameters 
     22   USE dom_oce        ! Ocean domain 
     23   USE sbc_oce        ! Surface boundary condition: ocean fields 
     24   USE sbc_ice        ! Surface boundary condition: ice fields 
    2925#if defined key_lim3 
    30    USE ice              ! LIM-3: ice variables 
    31    USE dom_ice          ! LIM-3: ice domain 
    32    USE limitd_me        ! LIM-3:  
     26   USE ice            ! LIM-3: ice variables 
     27   USE dom_ice        ! LIM-3: ice domain 
     28   USE limitd_me      ! LIM-3:  
    3329#else 
    34    USE ice_2            ! LIM2: ice variables 
    35    USE dom_ice_2        ! LIM2: ice domain 
     30   USE ice_2          ! LIM-2: ice variables 
     31   USE dom_ice_2      ! LIM-2: ice domain 
    3632#endif 
     33   USE lbclnk         ! Lateral Boundary Condition / MPP link 
     34   USE lib_mpp        ! MPP library 
     35   USE wrk_nemo       ! work arrays 
     36   USE in_out_manager ! I/O manager 
     37   USE prtctl         ! Print control 
     38   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3739 
    3840   IMPLICIT NONE 
    … …  
    4749#  include "vectopt_loop_substitute.h90" 
    4850   !!---------------------------------------------------------------------- 
    49    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     51   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    5052   !! $Id$ 
    5153   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    124126      REAL(wp) ::   zindb         ! ice (1) or not (0)       
    125127      REAL(wp) ::   zdummy        ! dummy argument 
     128      REAL(wp) ::   zintb, zintn  ! dummy argument 
    126129 
    127130      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength 
    … …  
    144147      REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12 
    145148      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity 
     149      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope: 
     150                                                              !   ocean surface (ssh_m) if ice is not embedded 
     151                                                              !   ice top surface if ice is embedded 
    146152       
    147153      !!------------------------------------------------------------------- 
    … …  
    150156      CALL wrk_alloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                ) 
    151157      CALL wrk_alloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds          ) 
    152       CALL wrk_alloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr                         ) 
     158      CALL wrk_alloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
    153159 
    154160#if  defined key_lim2 && ! defined key_lim2_vp 
    … …  
    231237      !  v_oce2: ocean v component on v points                         
    232238 
     239      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
     240          !                                             
     241          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
     242          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 
     243         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp      
     244          ! 
     245          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
     246          !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
     247         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
     248          ! 
     249         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0 
     250          ! 
     251      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
     252         zpice(:,:) = ssh_m(:,:) 
     253      ENDIF 
     254 
    233255      DO jj = k_j1+1, k_jpj-1 
    234256         DO ji = fs_2, fs_jpim1 
    … …  
    273295            ! include it later 
    274296 
    275             zdsshx =  ( ssh_m(ji+1,jj) - ssh_m(ji,jj) ) / e1u(ji,jj) 
    276             zdsshy =  ( ssh_m(ji,jj+1) - ssh_m(ji,jj) ) / e2v(ji,jj) 
     297            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj) 
     298            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj) 
    277299 
    278300            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx 
    … …  
    746768      CALL wrk_dealloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                ) 
    747769      CALL wrk_dealloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds          ) 
    748       CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr                         ) 
     770      CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
    749771 
    750772   END SUBROUTINE lim_rhg 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90

    r3294 r3625  
    1212   !!   'key_lim3' :                                   LIM sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    !!   lim_rst_opn     : open ice restart file 
    15    !!   lim_rst_write   : write of the restart file  
    16    !!   lim_rst_read    : read  the restart file  
    17    !!---------------------------------------------------------------------- 
    18    USE ice              ! sea-ice variables 
    19    USE par_ice          ! sea-ice parameters 
    20    USE dom_oce          ! ocean domain 
    21    USE sbc_oce          ! Surface boundary condition: ocean fields 
    22    USE sbc_ice          ! Surface boundary condition: ice fields 
    23    USE in_out_manager   ! I/O manager 
    24    USE iom              ! I/O library 
    25    USE lib_mpp          ! MPP library 
    26    USE wrk_nemo         ! work arrays 
     14   !!   lim_rst_opn   : open ice restart file 
     15   !!   lim_rst_write : write of the restart file  
     16   !!   lim_rst_read  : read  the restart file  
     17   !!---------------------------------------------------------------------- 
     18   USE ice            ! sea-ice variables 
     19   USE par_ice        ! sea-ice parameters 
     20   USE dom_oce        ! ocean domain 
     21   USE sbc_oce        ! Surface boundary condition: ocean fields 
     22   USE sbc_ice        ! Surface boundary condition: ice fields 
     23   USE in_out_manager ! I/O manager 
     24   USE iom            ! I/O library 
     25   USE lib_mpp        ! MPP library 
     26   USE wrk_nemo       ! work arrays 
     27   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2728 
    2829   IMPLICIT NONE 
    … …  
    3738 
    3839   !!---------------------------------------------------------------------- 
    39    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     40   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    4041   !! $Id$ 
    4142   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    402403                     zsmax = 4.5_wp 
    403404                     zsmin = 3.5_wp 
    404                      IF( sm_i(ji,jj,jl) .LT. zsmin ) THEN 
     405                     IF(     sm_i(ji,jj,jl) < zsmin ) THEN 
    405406                        zalpha = 1._wp 
    406                      ELSEIF( sm_i(ji,jj,jl) .LT.zsmax ) THEN 
     407                     ELSEIF( sm_i(ji,jj,jl) < zsmax ) THEN 
    407408                        zalpha = sm_i(ji,jj,jl) / ( zsmin - zsmax ) + zsmax / ( zsmax - zsmin ) 
    408409                     ELSE 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r3294 r3625  
    99   !!            3.3  ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea 
    1010   !!                 !                  + simplification of the ice-ocean stress calculation 
    11    !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
     11   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
     12   !!            3.5  ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass 
    1213   !!---------------------------------------------------------------------- 
    1314#if defined key_lim3 
    … …  
    3435   USE prtctl           ! Print control 
    3536   USE cpl_oasis3, ONLY : lk_cpl 
     37   USE oce,        ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n 
     38   USE dom_ice,    ONLY : tms 
     39   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3640 
    3741   IMPLICIT NONE 
    … …  
    4246   PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim 
    4347 
    44    REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice  
    4548   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values 
    4649   REAL(wp)  ::   rzero  = 0._wp     
    … …  
    5457#  include "vectopt_loop_substitute.h90" 
    5558   !!---------------------------------------------------------------------- 
    56    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     59   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    5760   !! $Id$ 
    5861   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    8689      !!              - qns     : sea heat flux: non solar 
    8790      !!              - emp     : freshwater budget: volume flux  
    88       !!              - emps    : freshwater budget: concentration/dillution  
     91      !!              - sfx     : salt flux  
    8992      !!              - fr_i    : ice fraction 
    9093      !!              - tn_ice  : sea-ice surface temperature 
    … …  
    97100      ! 
    98101      INTEGER  ::   ji, jj           ! dummy loop indices 
    99       INTEGER  ::   ierr             ! local integer 
    100       INTEGER  ::   ifvt, i1mfr, idfr               ! some switches 
    101       INTEGER  ::   iflt, ial, iadv, ifral, ifrdv 
    102       REAL(wp) ::   zinda, zfons, zpme              ! local scalars 
    103       REAL(wp), POINTER, DIMENSION(:,:) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes 
    104       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace 
     102      INTEGER  ::   ierr, ifvt, i1mfr, idfr           ! local integer 
     103      INTEGER  ::   iflt, ial , iadv , ifral, ifrdv   !   -      - 
     104      REAL(wp) ::   zinda, zemp, zemp_snow, zfmm      ! local scalars 
     105      REAL(wp) ::   zemp_snw                          !   -      - 
     106      REAL(wp) ::   zfcm1 , zfcm2                     !   -      - 
     107      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace 
    105108      !!--------------------------------------------------------------------- 
    106109       
    107       CALL wrk_alloc( jpi, jpj, zfcm1 , zfcm2 ) 
    108110      IF( lk_cpl )   CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp ) 
    109111 
    … …  
    139141 
    140142            !   computation the solar flux at ocean surface 
    141             zfcm1(ji,jj)   = pfrld(ji,jj) * qsr(ji,jj)  + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) 
     143            zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) 
    142144            ! fstric     Solar flux transmitted trough the ice 
    143145            ! qsr        Net short wave heat flux on free ocean 
    … …  
    146148 
    147149            !  computation the non solar heat flux at ocean surface 
    148             zfcm2(ji,jj) = - zfcm1(ji,jj)                  & 
    149                &           + iflt    * ( fscmbq(ji,jj) )   & ! total abl -> fscmbq is given to the ocean 
    150                ! fscmbq and ffltbif are obsolete 
    151                !              &           + iflt * ffltbif(ji,jj) !!! only if one category is used 
    152                &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   & 
    153                &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice   & 
    154                &           + fhmec(ji,jj)     & ! new contribution due to snow melt in ridging!! 
    155                &           + fheat_rpo(ji,jj) & ! contribution from ridge formation 
    156                &           + fheat_res(ji,jj) 
    157             ! fscmbq  Part of the solar radiation transmitted through the ice and going to the ocean 
    158             !         computed in limthd_zdf.F90 
    159             ! ffltbif Total heat content of the ice (brine pockets+ice) / delta_t 
     150            zfcm2 = - zfcm1                                                                     & ! ??? 
     151               &    + iflt    * fscmbq(ji,jj)                                                   & ! total ablation: heat given to the ocean 
     152               &    + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   & 
     153               &    + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice   & 
     154               &    + fhmec(ji,jj)                                                              & ! snow melt when ridging 
     155               &    + fheat_mec(ji,jj)                                                          & ! ridge formation 
     156               &    + fheat_res(ji,jj)                                                            ! residual heat flux 
    160157            ! qcmif   Energy needed to bring the ocean surface layer until its freezing (ok) 
    161158            ! qldif   heat balance of the lead (or of the open ocean) 
    162             ! qfvbq   i think this is wrong! 
    163             ! ---> Array used to store energy in case of total lateral ablation 
    164             ! qfvbq latent heat uptake/release after accretion/ablation 
    165             ! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead 
    166  
    167             IF ( num_sal == 2 ) zfcm2(ji,jj) = zfcm2(ji,jj) + & 
    168                fhbri(ji,jj) ! new contribution due to brine drainage  
    169  
    170             ! bottom radiative component is sent to the computation of the 
    171             ! oceanic heat flux 
    172             fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)      
     159            ! qfvbq   latent heat uptake/release after accretion/ablation 
     160            ! qdtcn   Energy from the turbulent oceanic heat flux heat flux coming in the lead 
     161 
     162            IF( num_sal == 2 )   zfcm2 = zfcm2 + fhbri(ji,jj)    ! add contribution due to brine drainage  
     163 
     164            ! bottom radiative component is sent to the computation of the oceanic heat flux 
     165            fsbbq(ji,jj) = ( 1._wp - ( ifvt + iflt ) ) * fscmbq(ji,jj)      
    173166 
    174167            ! used to compute the oceanic heat flux at the next time step 
    175             qsr(ji,jj) = zfcm1(ji,jj)                                       ! solar heat flux  
    176             qns(ji,jj) = zfcm2(ji,jj) - fdtcn(ji,jj)                        ! non solar heat flux 
     168            qsr(ji,jj) = zfcm1                                       ! solar heat flux  
     169            qns(ji,jj) = zfcm2 - fdtcn(ji,jj)                        ! non solar heat flux 
    177170            !                           ! fdtcn : turbulent oceanic heat flux 
    178171 
    179             !!gm   this IF prevents the vertorisation of the whole loop 
     172!!gm   this IF prevents the vertorisation of the whole loop 
    180173            IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 
    181174               WRITE(numout,*) ' lim_sbc : heat fluxes ' 
    182175               WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx) 
    183                WRITE(numout,*) ' zfcm1     : ', zfcm1(jiindx,jjindx) 
    184176               WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx) 
    185177               WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx) 
    186178               WRITE(numout,*) 
    187179               WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx) 
    188                WRITE(numout,*) ' zfcm2     : ', zfcm2(jiindx,jjindx) 
    189                WRITE(numout,*) ' zfcm1     : ', zfcm1(jiindx,jjindx) 
     180               WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
    190181               WRITE(numout,*) ' ifral     : ', ifral 
    191182               WRITE(numout,*) ' ial       : ', ial   
    … …  
    202193               WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
    203194               WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx) 
    204                WRITE(numout,*) ' fheat_rpo : ', fheat_rpo(jiindx,jjindx) 
     195               WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx) 
    205196               WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx) 
    206197               WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 
    207198            ENDIF 
    208             !!gm   end 
     199!!gm   end 
    209200         END DO 
    210201      END DO 
    … …  
    227218 
    228219            !  computing freshwater exchanges at the ice/ocean interface 
    229             zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction 
    230                &   + tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean 
    231                &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice 
    232                &   - rdmsnif(ji,jj) * r1_rdtice                       &   ! freshwaterflux due to snow melting  
    233                &   + fmmec(ji,jj)                                         ! snow falling when ridging 
    234  
    235  
    236             !  computing salt exchanges at the ice/ocean interface 
    237             !  sice should be the same as computed with the ice model 
    238             zfons =  ( soce_0(ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice  
    239             ! SOCE 
    240             zfons =  ( sss_m (ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice 
    241  
    242             !CT useless            !  salt flux for constant salinity 
    243             !CT useless            fsalt(ji,jj)      =  zfons / ( sss_m(ji,jj) + epsi16 ) + fsalt_res(ji,jj) 
    244             !  salt flux for variable salinity 
    245             zinda             = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 
    246             !  correcting brine and salt fluxes 
    247             fsbri(ji,jj)      =  zinda*fsbri(ji,jj) 
    248             !  converting the salt fluxes from ice to a freshwater flux from ocean 
    249             fsalt_res(ji,jj)  =  fsalt_res(ji,jj) / ( sss_m(ji,jj) + epsi16 ) 
    250             fseqv(ji,jj)      =  fseqv(ji,jj)     / ( sss_m(ji,jj) + epsi16 ) 
    251             fsbri(ji,jj)      =  fsbri(ji,jj)     / ( sss_m(ji,jj) + epsi16 ) 
    252             fsalt_rpo(ji,jj)  =  fsalt_rpo(ji,jj) / ( sss_m(ji,jj) + epsi16 ) 
    253  
    254             !  freshwater mass exchange (positive to the ice, negative for the ocean ?) 
    255             !  actually it's a salt flux (so it's minus freshwater flux) 
    256             !  if sea ice grows, zfons is positive, fsalt also 
    257             !  POSITIVE SALT FLUX FROM THE ICE TO THE OCEAN 
    258             !  POSITIVE FRESHWATER FLUX FROM THE OCEAN TO THE ICE [kg.m-2.s-1] 
    259  
    260             emp(ji,jj) = - zpme  
     220            zemp =   emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction 
     221               &   - tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean 
     222               &   + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice 
     223               &   - fmmec(ji,jj)                                         ! snow falling when ridging 
     224 
     225            ! mass flux at the ocean/ice interface (sea ice fraction) 
     226            zemp_snw = rdm_snw(ji,jj) * r1_rdtice                         ! snow melting = pure water that enters the ocean 
     227            zfmm     = rdm_ice(ji,jj) * r1_rdtice                         ! Freezing minus mesting   
     228 
     229            emp(ji,jj) = zemp + zemp_snw + zfmm  ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
     230             
     231            !  correcting brine salt fluxes   (zinda = 1  if pfrld=1 , =0 otherwise) 
     232            zinda        = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 
     233            sfx_bri(ji,jj) = zinda * sfx_bri(ji,jj) 
    261234         END DO 
    262235      END DO 
    263236 
     237      !------------------------------------------! 
     238      !      salt flux at the ocean surface      ! 
     239      !------------------------------------------! 
     240 
    264241      IF( num_sal == 2 ) THEN      ! variable ice salinity: brine drainage included in the salt flux 
    265          emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:) 
     242         sfx(:,:) = sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) + sfx_bri(:,:) 
    266243      ELSE                         ! constant ice salinity: 
    267          emps(:,:) =              fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:) 
     244         sfx(:,:) = sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) 
     245      ENDIF 
     246      !-----------------------------------------------! 
     247      !   mass of snow and ice per unit area          ! 
     248      !-----------------------------------------------! 
     249      IF( nn_ice_embd /= 0 ) THEN                               ! embedded sea-ice (mass required) 
     250         snwice_mass_b(:,:) = snwice_mass(:,:)                  ! save mass from the previous ice time step 
     251         !                                                      ! new mass per unit area 
     252         snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )  
     253         !                                                      ! time evolution of snow+ice mass 
     254         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice 
    268255      ENDIF 
    269256 
    … …  
    285272      IF(ln_ctl) THEN 
    286273         CALL prt_ctl( tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns , clinfo2=' qns     : ' ) 
    287          CALL prt_ctl( tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps, clinfo2=' emps    : ' ) 
     274         CALL prt_ctl( tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=sfx , clinfo2=' sfx     : ' ) 
    288275         CALL prt_ctl( tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ' ) 
    289276         CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl ) 
    290277      ENDIF 
    291278      ! 
    292       CALL wrk_dealloc( jpi, jpj, zfcm1 , zfcm2 ) 
    293279      IF( lk_cpl )   CALL wrk_dealloc( jpi, jpj, jpl, zalb, zalbp ) 
    294280      !  
    … …  
    383369      !!------------------------------------------------------------------- 
    384370      ! 
     371      INTEGER  ::   ji, jj                          ! dummy loop indices 
     372      REAL(wp) ::   zcoefu, zcoefv, zcoeff          ! local scalar 
    385373      IF(lwp) WRITE(numout,*) 
    386374      IF(lwp) WRITE(numout,*) 'lim_sbc_init : LIM-3 sea-ice - surface boundary condition' 
    … …  
    389377      !                                      ! allocate lim_sbc array 
    390378      IF( lim_sbc_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' ) 
    391       ! 
    392       r1_rdtice = 1. / rdt_ice 
    393379      ! 
    394380      soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case 
    … …  
    402388         END WHERE 
    403389      ENDIF 
     390      !                                      ! embedded sea ice 
     391      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass 
     392         snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
     393         snwice_mass_b(:,:) = snwice_mass(:,:) 
     394      ELSE 
     395         snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges 
     396         snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges 
     397      ENDIF 
     398      IF( nn_ice_embd == 2  .AND.         &  ! full embedment (case 2) & no restart 
     399         &  .NOT. ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area 
     400         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
     401         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 
     402         ! 
     403         ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n]  
     404         !       which were previously set in domvvl 
     405         IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only??? 
     406            DO jj = 1, jpjm1 
     407               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible 
     408                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 
     409                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 
     410                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 
     411                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     & 
     412                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 
     413                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     & 
     414                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 
     415                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     & 
     416                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 
     417                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     & 
     418                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 
     419               END DO 
     420            END DO 
     421            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. ) 
     422            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. ) 
     423            DO jj = 1, jpjm1 
     424               DO ji = 1, jpim1      ! NO Vector Opt. 
     425                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   & 
     426                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     & 
     427                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     & 
     428                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 
     429               END DO 
     430            END DO 
     431            CALL lbc_lnk( sshf_n, 'F', 1. ) 
     432          ENDIF 
     433      ENDIF 
    404434      ! 
    405435   END SUBROUTINE lim_sbc_init 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90

    r2715 r3625  
    22   !!====================================================================== 
    33   !!                       ***  MODULE limtab   *** 
    4    !!   LIM : transform 1D (2D) array to a 2D (1D) table 
     4   !!   LIM ice model : transform 1D (2D) array to a 2D (1D) table 
    55   !!====================================================================== 
    66#if defined key_lim3 
    … …  
    2020 
    2121   !!---------------------------------------------------------------------- 
    22    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
     22   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2010) 
    2323   !! $Id$ 
    2424   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r3294 r3625  
    88   !!            3.0  ! 2005-11 (M. Vancoppenolle)  LIM-3 : Multi-layer thermodynamics + salinity variations 
    99   !!             -   ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec, lim_thd_con_dh and lim_thd_con_dif 
    10    !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif 
     10   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw 
    1111   !!            3.3  ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas) 
    1212   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
    … …  
    1616   !!   'key_lim3'                                      LIM3 sea-ice model 
    1717   !!---------------------------------------------------------------------- 
    18    !!   lim_thd        : thermodynamic of sea ice 
    19    !!   lim_thd_init   : initialisation of sea-ice thermodynamic 
     18   !!   lim_thd       : thermodynamic of sea ice 
     19   !!   lim_thd_init  : initialisation of sea-ice thermodynamic 
    2020   !!---------------------------------------------------------------------- 
    21    USE phycst          ! physical constants 
    22    USE dom_oce         ! ocean space and time domain variables 
    23    USE ice             ! LIM: sea-ice variables 
    24    USE par_ice         ! LIM: sea-ice parameters 
    25    USE sbc_oce         ! Surface boundary condition: ocean fields 
    26    USE sbc_ice         ! Surface boundary condition: ice fields 
    27    USE thd_ice         ! LIM thermodynamic sea-ice variables 
    28    USE dom_ice         ! LIM sea-ice domain 
    29    USE domvvl          ! domain: variable volume level 
    30    USE limthd_dif      ! LIM: thermodynamics, vertical diffusion 
    31    USE limthd_dh       ! LIM: thermodynamics, ice and snow thickness variation 
    32    USE limthd_sal      ! LIM: thermodynamics, ice salinity 
    33    USE limthd_ent      ! LIM: thermodynamics, ice enthalpy redistribution 
    34    USE limtab          ! LIM: 1D <==> 2D transformation 
    35    USE limvar          ! LIM: sea-ice variables 
    36    USE lbclnk          ! lateral boundary condition - MPP links 
    37    USE lib_mpp         ! MPP library 
    38    USE wrk_nemo        ! work arrays 
    39    USE in_out_manager  ! I/O manager 
    40    USE prtctl          ! Print control 
     21   USE phycst         ! physical constants 
     22   USE dom_oce        ! ocean space and time domain variables 
     23   USE ice            ! LIM: sea-ice variables 
     24   USE par_ice        ! LIM: sea-ice parameters 
     25   USE sbc_oce        ! Surface boundary condition: ocean fields 
     26   USE sbc_ice        ! Surface boundary condition: ice fields 
     27   USE thd_ice        ! LIM thermodynamic sea-ice variables 
     28   USE dom_ice        ! LIM sea-ice domain 
     29   USE domvvl         ! domain: variable volume level 
     30   USE limthd_dif     ! LIM: thermodynamics, vertical diffusion 
     31   USE limthd_dh      ! LIM: thermodynamics, ice and snow thickness variation 
     32   USE limthd_sal     ! LIM: thermodynamics, ice salinity 
     33   USE limthd_ent     ! LIM: thermodynamics, ice enthalpy redistribution 
     34   USE limtab         ! LIM: 1D <==> 2D transformation 
     35   USE limvar         ! LIM: sea-ice variables 
     36   USE lbclnk         ! lateral boundary condition - MPP links 
     37   USE lib_mpp        ! MPP library 
     38   USE wrk_nemo       ! work arrays 
     39   USE in_out_manager ! I/O manager 
     40   USE prtctl         ! Print control 
     41   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    4142 
    4243   IMPLICIT NONE 
    … …  
    110111                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * nlay_i 
    111112                  !0 if no ice and 1 if yes 
    112                   zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i(ji,jj,jl) ) )  
     113                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) ) )  
    113114                  !convert units ! very important that this line is here 
    114115                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb  
    … …  
    122123                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * nlay_s 
    123124                  !0 if no ice and 1 if yes 
    124                   zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s(ji,jj,jl) ) )  
     125                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) ) )  
    125126                  !convert units ! very important that this line is here 
    126127                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb  
    … …  
    140141      ffltbif(:,:) = 0.e0   ! linked with fstric 
    141142      qfvbq  (:,:) = 0.e0   ! linked with fstric 
    142       rdmsnif(:,:) = 0.e0   ! variation of snow mass per unit area 
    143       rdmicif(:,:) = 0.e0   ! variation of ice mass per unit area 
     143      rdm_snw(:,:) = 0.e0   ! variation of snow mass per unit area 
     144      rdm_ice(:,:) = 0.e0   ! variation of ice mass per unit area 
    144145      hicifp (:,:) = 0.e0   ! daily thermodynamic ice production.  
    145       fsbri  (:,:) = 0.e0   ! brine flux contribution to salt flux to the ocean 
     146      sfx_bri(:,:) = 0.e0   ! brine flux contribution to salt flux to the ocean 
    146147      fhbri  (:,:) = 0.e0   ! brine flux contribution to heat flux to the ocean 
    147       fseqv  (:,:) = 0.e0   ! equivalent salt flux to the ocean due to ice/growth decay 
     148      sfx_thd(:,:) = 0.e0   ! equivalent salt flux to the ocean due to ice/growth decay 
    148149 
    149150      !----------------------------------- 
    … …  
    273274            CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb), fr2_i0          , jpi, jpj, npb(1:nbpb) ) 
    274275            CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
    275  
    276276#if ! defined key_coupled 
    277             CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) 
    278             CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl)   , jpi, jpj, npb(1:nbpb) ) 
     277            CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) ) 
     278            CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
    279279#endif 
    280  
    281             CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl)   , jpi, jpj, npb(1:nbpb) ) 
    282             CALL tab_2d_1d( nbpb, t_bo_b     (1:nbpb), t_bo       , jpi, jpj, npb(1:nbpb) ) 
    283             CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip    , jpi, jpj, npb(1:nbpb) )  
    284             CALL tab_2d_1d( nbpb, fbif_1d    (1:nbpb), fbif       , jpi, jpj, npb(1:nbpb) ) 
    285             CALL tab_2d_1d( nbpb, qldif_1d   (1:nbpb), qldif      , jpi, jpj, npb(1:nbpb) ) 
    286             CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb), rdmicif    , jpi, jpj, npb(1:nbpb) ) 
    287             CALL tab_2d_1d( nbpb, rdmsnif_1d (1:nbpb), rdmsnif    , jpi, jpj, npb(1:nbpb) ) 
    288             CALL tab_2d_1d( nbpb, dmgwi_1d   (1:nbpb), dmgwi      , jpi, jpj, npb(1:nbpb) ) 
    289             CALL tab_2d_1d( nbpb, qlbbq_1d   (1:nbpb), zqlbsbq    , jpi, jpj, npb(1:nbpb) ) 
    290  
    291             CALL tab_2d_1d( nbpb, fseqv_1d   (1:nbpb), fseqv      , jpi, jpj, npb(1:nbpb) ) 
    292             CALL tab_2d_1d( nbpb, fsbri_1d   (1:nbpb), fsbri      , jpi, jpj, npb(1:nbpb) ) 
    293             CALL tab_2d_1d( nbpb, fhbri_1d   (1:nbpb), fhbri      , jpi, jpj, npb(1:nbpb) ) 
    294             CALL tab_2d_1d( nbpb, fstbif_1d  (1:nbpb), fstric     , jpi, jpj, npb(1:nbpb) ) 
    295             CALL tab_2d_1d( nbpb, qfvbq_1d   (1:nbpb), qfvbq      , jpi, jpj, npb(1:nbpb) ) 
     280            CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) ) 
     281            CALL tab_2d_1d( nbpb, t_bo_b     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) ) 
     282            CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip         , jpi, jpj, npb(1:nbpb) )  
     283            CALL tab_2d_1d( nbpb, fbif_1d    (1:nbpb), fbif            , jpi, jpj, npb(1:nbpb) ) 
     284            CALL tab_2d_1d( nbpb, qldif_1d   (1:nbpb), qldif           , jpi, jpj, npb(1:nbpb) ) 
     285            CALL tab_2d_1d( nbpb, rdm_ice_1d (1:nbpb), rdm_ice         , jpi, jpj, npb(1:nbpb) ) 
     286            CALL tab_2d_1d( nbpb, rdm_snw_1d (1:nbpb), rdm_snw         , jpi, jpj, npb(1:nbpb) ) 
     287            CALL tab_2d_1d( nbpb, dmgwi_1d   (1:nbpb), dmgwi           , jpi, jpj, npb(1:nbpb) ) 
     288            CALL tab_2d_1d( nbpb, qlbbq_1d   (1:nbpb), zqlbsbq         , jpi, jpj, npb(1:nbpb) ) 
     289 
     290            CALL tab_2d_1d( nbpb, sfx_thd_1d (1:nbpb), sfx_thd         , jpi, jpj, npb(1:nbpb) ) 
     291            CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) ) 
     292            CALL tab_2d_1d( nbpb, fhbri_1d   (1:nbpb), fhbri           , jpi, jpj, npb(1:nbpb) ) 
     293            CALL tab_2d_1d( nbpb, fstbif_1d  (1:nbpb), fstric          , jpi, jpj, npb(1:nbpb) ) 
     294            CALL tab_2d_1d( nbpb, qfvbq_1d   (1:nbpb), qfvbq           , jpi, jpj, npb(1:nbpb) ) 
    296295 
    297296            !-------------------------------- 
    … …  
    331330            !-------------------------------- 
    332331 
    333             CALL tab_1d_2d( nbpb, at_i        , npb, at_i_b(1:nbpb), jpi, jpj ) 
    334             CALL tab_1d_2d( nbpb, ht_i(:,:,jl), npb, ht_i_b(1:nbpb), jpi, jpj ) 
    335             CALL tab_1d_2d( nbpb, ht_s(:,:,jl), npb, ht_s_b(1:nbpb), jpi, jpj ) 
    336             CALL tab_1d_2d( nbpb, a_i (:,:,jl), npb, a_i_b(1:nbpb) , jpi, jpj ) 
    337             CALL tab_1d_2d( nbpb, t_su(:,:,jl), npb, t_su_b(1:nbpb), jpi, jpj ) 
    338             CALL tab_1d_2d( nbpb, sm_i(:,:,jl), npb, sm_i_b(1:nbpb), jpi, jpj ) 
    339  
     332               CALL tab_1d_2d( nbpb, at_i          , npb, at_i_b    (1:nbpb)   , jpi, jpj ) 
     333               CALL tab_1d_2d( nbpb, ht_i(:,:,jl)  , npb, ht_i_b    (1:nbpb)   , jpi, jpj ) 
     334               CALL tab_1d_2d( nbpb, ht_s(:,:,jl)  , npb, ht_s_b    (1:nbpb)   , jpi, jpj ) 
     335               CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_b     (1:nbpb)   , jpi, jpj ) 
     336               CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_b    (1:nbpb)   , jpi, jpj ) 
     337               CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_b    (1:nbpb)   , jpi, jpj ) 
    340338            DO jk = 1, nlay_s 
    341                CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b(1:nbpb,jk), jpi, jpj) 
    342                CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_b(1:nbpb,jk), jpi, jpj) 
     339               CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b     (1:nbpb,jk), jpi, jpj) 
     340               CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_b     (1:nbpb,jk), jpi, jpj) 
    343341            END DO 
    344  
    345342            DO jk = 1, nlay_i 
    346                CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b(1:nbpb,jk), jpi, jpj) 
    347                CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b(1:nbpb,jk), jpi, jpj) 
    348                CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b(1:nbpb,jk), jpi, jpj) 
     343               CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b     (1:nbpb,jk), jpi, jpj) 
     344               CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b     (1:nbpb,jk), jpi, jpj) 
     345               CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b     (1:nbpb,jk), jpi, jpj) 
    349346            END DO 
    350  
    351             CALL tab_1d_2d( nbpb, fstric , npb, fstbif_1d (1:nbpb), jpi, jpj ) 
    352             CALL tab_1d_2d( nbpb, qldif  , npb, qldif_1d  (1:nbpb), jpi, jpj ) 
    353             CALL tab_1d_2d( nbpb, qfvbq  , npb, qfvbq_1d  (1:nbpb), jpi, jpj ) 
    354             CALL tab_1d_2d( nbpb, rdmicif, npb, rdmicif_1d(1:nbpb), jpi, jpj ) 
    355             CALL tab_1d_2d( nbpb, rdmsnif, npb, rdmsnif_1d(1:nbpb), jpi, jpj ) 
    356             CALL tab_1d_2d( nbpb, dmgwi  , npb, dmgwi_1d  (1:nbpb), jpi, jpj ) 
    357             CALL tab_1d_2d( nbpb, rdvosif, npb, dvsbq_1d  (1:nbpb), jpi, jpj ) 
    358             CALL tab_1d_2d( nbpb, rdvobif, npb, dvbbq_1d  (1:nbpb), jpi, jpj ) 
    359             CALL tab_1d_2d( nbpb, fdvolif, npb, dvlbq_1d  (1:nbpb), jpi, jpj ) 
    360             CALL tab_1d_2d( nbpb, rdvonif, npb, dvnbq_1d  (1:nbpb), jpi, jpj )  
    361             CALL tab_1d_2d( nbpb, fseqv  , npb, fseqv_1d  (1:nbpb), jpi, jpj ) 
     347               CALL tab_1d_2d( nbpb, fstric        , npb, fstbif_1d (1:nbpb)   , jpi, jpj ) 
     348               CALL tab_1d_2d( nbpb, qldif         , npb, qldif_1d  (1:nbpb)   , jpi, jpj ) 
     349               CALL tab_1d_2d( nbpb, qfvbq         , npb, qfvbq_1d  (1:nbpb)   , jpi, jpj ) 
     350               CALL tab_1d_2d( nbpb, rdm_ice       , npb, rdm_ice_1d(1:nbpb)   , jpi, jpj ) 
     351               CALL tab_1d_2d( nbpb, rdm_snw       , npb, rdm_snw_1d(1:nbpb)   , jpi, jpj ) 
     352               CALL tab_1d_2d( nbpb, dmgwi         , npb, dmgwi_1d  (1:nbpb)   , jpi, jpj ) 
     353               CALL tab_1d_2d( nbpb, rdvosif       , npb, dvsbq_1d  (1:nbpb)   , jpi, jpj ) 
     354               CALL tab_1d_2d( nbpb, rdvobif       , npb, dvbbq_1d  (1:nbpb)   , jpi, jpj ) 
     355               CALL tab_1d_2d( nbpb, fdvolif       , npb, dvlbq_1d  (1:nbpb)   , jpi, jpj ) 
     356               CALL tab_1d_2d( nbpb, rdvonif       , npb, dvnbq_1d  (1:nbpb)   , jpi, jpj )  
     357               CALL tab_1d_2d( nbpb, sfx_thd       , npb, sfx_thd_1d(1:nbpb)   , jpi, jpj ) 
    362358            ! 
    363359            IF( num_sal == 2 ) THEN 
    364                CALL tab_1d_2d( nbpb, fsbri, npb, fsbri_1d(1:nbpb), jpi, jpj ) 
    365                CALL tab_1d_2d( nbpb, fhbri, npb, fhbri_1d(1:nbpb), jpi, jpj ) 
     360               CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj ) 
     361               CALL tab_1d_2d( nbpb, fhbri         , npb, fhbri_1d  (1:nbpb)   , jpi, jpj ) 
    366362            ENDIF 
    367363            ! 
    368             !+++++ 
    369             !temporary stuff for a dummy version 
     364            !+++++       temporary stuff for a dummy version 
    370365            CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb)      , jpi, jpj ) 
    371366            CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb)      , jpi, jpj ) 
    … …  
    389384      ! 5.1) Ice heat content               
    390385      !------------------------ 
    391       ! Enthalpies are global variables we have to readjust the units 
     386      ! Enthalpies are global variables we have to readjust the units (heat content in 10^9 Joules) 
    392387      zcoef = 1._wp / ( unit_fac * REAL( nlay_i ) ) 
    393388      DO jl = 1, jpl 
    394389         DO jk = 1, nlay_i 
    395             ! Multiply by volume, divide by nlayers so that heat content in 10^9 Joules 
    396390            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) * zcoef 
    397391         END DO 
    … …  
    401395      ! 5.2) Snow heat content               
    402396      !------------------------ 
    403       ! Enthalpies are global variables we have to readjust the units 
     397      ! Enthalpies are global variables we have to readjust the units (heat content in 10^9 Joules) 
    404398      zcoef = 1._wp / ( unit_fac * REAL( nlay_s ) ) 
    405399      DO jl = 1, jpl 
    406400         DO jk = 1, nlay_s 
    407             ! Multiply by volume, so that heat content in 10^9 Joules 
    408401            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) * zcoef 
    409402         END DO 
    … …  
    419412      !-------------------------------------------- 
    420413      d_v_i_thd(:,:,:) = v_i      (:,:,:) - old_v_i(:,:,:)    ! ice volumes  
    421       dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) / rdt_ice * 86400.0 
     414      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    422415 
    423416      IF( con_i )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 
    … …  
    488481      ! 
    489482      IF(lwp) WRITE(numout,*) ' lim_thd_glohec ' 
    490       IF(lwp) WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) / rdt_ice 
    491       IF(lwp) WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) / rdt_ice 
    492       IF(lwp) WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) / rdt_ice 
     483      IF(lwp) WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice 
     484      IF(lwp) WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice 
     485      IF(lwp) WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice 
    493486      ! 
    494487   END SUBROUTINE lim_thd_glohec 
    … …  
    538531      !-------------------- 
    539532      DO ji = kideb, kiut 
    540          cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 
     533         cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 
    541534      END DO 
    542535 
    … …  
    597590            WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 
    598591            WRITE(numout,*) ' surf_error : ', surf_error(ji,jl) 
    599             WRITE(numout,*) ' dq_i       : ', - dq_i(ji,jl) / rdt_ice 
     592            WRITE(numout,*) ' dq_i       : ', - dq_i(ji,jl) * r1_rdtice 
    600593            WRITE(numout,*) ' Fdt        : ', sum_fluxq(ji,jl) 
    601594            WRITE(numout,*) 
    … …  
    631624            WRITE(numout,*) 
    632625            WRITE(numout,*) ' Layer by layer ... ' 
    633             WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice 
     626            WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) * r1_rdtice 
    634627            WRITE(numout,*) ' dfc_snow  : ', fc_s(ji,1) - fc_s(ji,0) 
    635628            DO jk = 1, nlay_i 
    636629               WRITE(numout,*) ' layer  : ', jk 
    637                WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice   
     630               WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) * r1_rdtice   
    638631               WRITE(numout,*) ' radab  : ', radab(ji,jk) 
    639632               WRITE(numout,*) ' dfc_i  : ', fc_i(ji,jk) - fc_i(ji,jk-1) 
    … …  
    681674         fatm      (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji)                       ! total heat flux 
    682675         sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(zji,zjj,jl)  
    683          cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 
     676         cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 
    684677      END DO 
    685678 
    … …  
    688681      !-------------------- 
    689682      DO ji = kideb, kiut 
    690          cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 
     683         cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 
    691684      END DO 
    692685 
    … …  
    722715            WRITE(numout,*) ' * ' 
    723716            WRITE(numout,*) ' Ftotal     : ', sum_fluxq(ji,jl) 
    724             WRITE(numout,*) ' dq_t       : ', - dq_i(ji,jl) / rdt_ice 
    725             WRITE(numout,*) ' dq_i       : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) / rdt_ice 
    726             WRITE(numout,*) ' dq_s       : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice 
     717            WRITE(numout,*) ' dq_t       : ', - dq_i(ji,jl) * r1_rdtice 
     718            WRITE(numout,*) ' dq_i       : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) * r1_rdtice 
     719            WRITE(numout,*) ' dq_s       : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) * r1_rdtice 
    727720            WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 
    728721            WRITE(numout,*) ' * ' 
    … …  
    734727            WRITE(numout,*) ' * ' 
    735728            WRITE(numout,*) ' Heat contents --- : ' 
    736             WRITE(numout,*) ' qt_s_in    : ', qt_s_in(ji,jl) / rdt_ice 
    737             WRITE(numout,*) ' qt_i_in    : ', qt_i_in(ji,jl) / rdt_ice 
    738             WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) / rdt_ice 
    739             WRITE(numout,*) ' qt_s_fin   : ', qt_s_fin(ji,jl) / rdt_ice 
    740             WRITE(numout,*) ' qt_i_fin   : ', qt_i_fin(ji,jl) / rdt_ice 
    741             WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) / rdt_ice 
     729            WRITE(numout,*) ' qt_s_in    : ', qt_s_in(ji,jl) * r1_rdtice 
     730            WRITE(numout,*) ' qt_i_in    : ', qt_i_in(ji,jl) * r1_rdtice 
     731            WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) * r1_rdtice 
     732            WRITE(numout,*) ' qt_s_fin   : ', qt_s_fin(ji,jl) * r1_rdtice 
     733            WRITE(numout,*) ' qt_i_fin   : ', qt_i_fin(ji,jl) * r1_rdtice 
     734            WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) * r1_rdtice 
    742735            WRITE(numout,*) ' * ' 
    743736            WRITE(numout,*) ' Ice variables --- : ' 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r3294 r3625  
    77   !!                 ! 2005-06 (M. Vancoppenolle) 3D version  
    88   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif & rdmicif 
    9    !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
     9   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
     10   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes  
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_lim3 
    … …  
    1314   !!   'key_lim3'                                      LIM3 sea-ice model 
    1415   !!---------------------------------------------------------------------- 
    15    !!   lim_thd_dh  : vertical accr./abl. and lateral ablation of sea ice 
    16    !!---------------------------------------------------------------------- 
    17    USE par_oce          ! ocean parameters 
    18    USE phycst           ! physical constants (OCE directory)  
    19    USE sbc_oce          ! Surface boundary condition: ocean fields 
    20    USE ice              ! LIM variables 
    21    USE par_ice          ! LIM parameters 
    22    USE thd_ice          ! LIM thermodynamics 
    23    USE in_out_manager   ! I/O manager 
    24    USE lib_mpp          ! MPP library 
    25    USE wrk_nemo         ! work arrays 
     16   !!   lim_thd_dh    : vertical accr./abl. and lateral ablation of sea ice 
     17   !!---------------------------------------------------------------------- 
     18   USE par_oce        ! ocean parameters 
     19   USE phycst         ! physical constants (OCE directory)  
     20   USE sbc_oce        ! Surface boundary condition: ocean fields 
     21   USE ice            ! LIM variables 
     22   USE par_ice        ! LIM parameters 
     23   USE thd_ice        ! LIM thermodynamics 
     24   USE in_out_manager ! I/O manager 
     25   USE lib_mpp        ! MPP library 
     26   USE wrk_nemo       ! work arrays 
     27   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2628 
    2729   IMPLICIT NONE 
    … …  
    3739 
    3840   !!---------------------------------------------------------------------- 
    39    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
     41   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    4042   !! $Id$ 
    4143   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    7173      !!  
    7274      INTEGER  ::   ji , jk        ! dummy loop indices 
    73       INTEGER  ::   zji, zjj       ! 2D corresponding indices to ji 
     75      INTEGER  ::   ii, ij       ! 2D corresponding indices to ji 
    7476      INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow 
    7577      INTEGER  ::   isnowic        ! snow ice formation not 
    … …  
    102104      REAL(wp), POINTER, DIMENSION(:) ::   zfmass_i    !  
    103105 
    104       REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel     ! snow melt  
    105       REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre     ! snow precipitation  
    106       REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub     ! snow sublimation 
    107       REAL(wp), POINTER, DIMENSION(:) ::   zfsalt_melt   ! salt flux due to ice melt 
     106      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt  
     107      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre   ! snow precipitation  
     108      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub   ! snow sublimation 
     109      REAL(wp), POINTER, DIMENSION(:) ::   zsfx_melt   ! salt flux due to ice melt 
    108110 
    109111      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah 
    … …  
    126128 
    127129      CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 
    128       CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfsalt_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
     130      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zsfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
    129131      CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i ) 
    130132      CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) 
    131133 
    132       zfsalt_melt(:) = 0._wp 
    133       ftotal_fin(:)   = 0._wp 
    134       zfdt_init(:)    = 0._wp 
    135       zfdt_final(:)   = 0._wp 
     134      zsfx_melt (:) = 0._wp 
     135      ftotal_fin(:) = 0._wp 
     136      zfdt_init (:) = 0._wp 
     137      zfdt_final(:) = 0._wp 
    136138 
    137139      DO ji = kideb, kiut 
    … …  
    145147      ! 
    146148      DO ji = kideb, kiut 
    147          isnow         = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) ) 
    148          ztfs(ji)      = isnow * rtt + ( 1.0 - isnow ) * rtt 
    149          z_f_surf(ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
    150          z_f_surf(ji)  = MAX( zzero , z_f_surf(ji) ) * MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) 
     149         isnow         = INT(  1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s_b(ji) )  ) ) 
     150         ztfs     (ji) = isnow * rtt + ( 1.0 - isnow ) * rtt 
     151         z_f_surf (ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
     152         z_f_surf (ji) = MAX(  zzero , z_f_surf(ji)  ) * MAX(  zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) 
    151153         zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice 
    152154      END DO ! ji 
    … …  
    240242         zhsnew         =  ht_s_b(ji) + dh_s_tot(ji) 
    241243         ! If snow is still present zhn = 1, else zhn = 0 
    242          zhn            =  1.0 - MAX( zzero , SIGN( zone , - zhsnew ) ) 
     244         zhn            =  1.0 - MAX(  zzero , SIGN( zone , - zhsnew ) ) 
    243245         ht_s_b(ji)     =  MAX( zzero , zhsnew ) 
    244246         ! Volume and mass variations of snow 
    245          dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) ) 
     247         dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) 
    246248         dvsbq_1d  (ji) =  MIN( zzero, dvsbq_1d(ji) ) 
    247          rdmsnif_1d(ji) =  rdmsnif_1d(ji) + rhosn * dvsbq_1d(ji) 
     249         rdm_snw_1d(ji) =  rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 
    248250      END DO ! ji 
    249251 
    … …  
    253255      DO ji = kideb, kiut  
    254256         dh_i_surf(ji) =  0._wp 
    255          z_f_surf (ji) =  zqfont_su(ji) / rdt_ice ! heat conservation test 
     257         z_f_surf (ji) =  zqfont_su(ji) * r1_rdtice  ! heat conservation test 
    256258         zdq_i    (ji) =  0._wp 
    257259      END DO ! ji 
    … …  
    262264            zdeltah  (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 
    263265            !                                                    ! recompute heat available 
    264             zqfont_su(ji)    = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)  
     266            zqfont_su(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)  
    265267            !                                                    ! melt of layer jk cannot be higher than its thickness 
    266268            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 
    267269            !                                                    ! update surface melt 
    268             dh_i_surf(ji)    = dh_i_surf(ji) + zdeltah(ji,jk)  
     270            dh_i_surf(ji   ) = dh_i_surf(ji) + zdeltah(ji,jk)  
    269271            !                                                    ! for energy conservation 
    270             zdq_i    (ji)    = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 
     272            zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    271273            ! 
    272             ! contribution to ice-ocean salt flux  
    273             zji = MOD( npb(ji) - 1 , jpi ) + 1 
    274             zjj =    ( npb(ji) - 1 ) / jpi + 1 
    275             zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji) ) * a_i_b(ji)    & 
    276                &                              * MIN( zdeltah(ji,jk) , 0.e0 ) * rhoic / rdt_ice  
     274            !                                                    ! contribution to ice-ocean salt flux  
     275            zsfx_melt(ji)  = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice  
    277276         END DO 
    278277      END DO 
    … …  
    290289            IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN! 
    291290               WRITE(numout,*) ' ALERTE heat loss for surface melt ' 
    292                WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl 
     291               WRITE(numout,*) ' ii, ij, jl :', ii, ij, jl 
    293292               WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji) 
    294293               WRITE(numout,*) ' z_f_surf     : ', z_f_surf(ji) 
    … …  
    299298               WRITE(numout,*) ' qlbbq_1d     : ', qlbbq_1d(ji) 
    300299               WRITE(numout,*) ' s_i_new      : ', s_i_new(ji) 
    301                WRITE(numout,*) ' sss_m        : ', sss_m(zji,zjj) 
     300               WRITE(numout,*) ' sss_m        : ', sss_m(ii,ij) 
    302301            ENDIF 
    303302         END DO 
    … …  
    338337         DO ji = kideb, kiut 
    339338            ! In case of disparition of the snow, we have to update the snow temperatures 
    340             zhisn  =  MAX( zzero , SIGN( zone, - ht_s_b(ji) ) ) 
     339            zhisn  =  MAX(  zzero , SIGN( zone, - ht_s_b(ji) ) ) 
    341340            t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt 
    342341            q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk) 
    … …  
    358357      ! 4.1 Basal growth - (a) salinity not varying in time  
    359358      !----------------------------------------------------- 
    360       IF(  num_sal /= 2  .AND.  num_sal /= 4  ) THEN 
    361          DO ji = kideb, kiut 
    362             IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0.0  ) THEN 
     359      IF(  num_sal /= 2  ) THEN   ! ice salinity constant in time 
     360         DO ji = kideb, kiut 
     361            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp  ) THEN 
    363362               s_i_new(ji)         =  sm_i_b(ji) 
    364363               ! Melting point in K 
    … …  
    371370                  &                           - rcp  * ( ztmelts - rtt )                                 ) 
    372371               ! Basal growth rate = - F*dt / q 
    373                dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
     372               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
    374373            ENDIF 
    375374         END DO 
    … …  
    379378      ! 4.1 Basal growth - (b) salinity varying in time  
    380379      !------------------------------------------------- 
    381       IF(  num_sal == 2 .OR.  num_sal == 4  ) THEN 
    382          ! the growth rate (dh_i_bott) is function of the new ice 
    383          ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice  
    384          ! salinity (snewice). snewice depends on dh_i_bott 
    385          ! it converges quickly, so, no problem 
     380      IF(  num_sal == 2  ) THEN 
     381         ! the growth rate (dh_i_bott) is function of the new ice heat content (q_i_b(nlay_i+1)).  
     382         ! q_i_b depends on the new ice salinity (snewice).  
     383         ! snewice depends on dh_i_bott ; it converges quickly, so, no problem 
    386384         ! See Vancoppenolle et al., OM08 for more info on this 
    387385 
    … …  
    394392            DO ji = kideb, kiut 
    395393               IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN 
    396                   zji = MOD( npb(ji) - 1, jpi ) + 1 
    397                   zjj = ( npb(ji) - 1 ) / jpi + 1 
     394                  ii = MOD( npb(ji) - 1, jpi ) + 1 
     395                  ij = ( npb(ji) - 1 ) / jpi + 1 
    398396                  ! Melting point in K 
    399397                  ztmelts             =   - tmut * s_i_new(ji) + rtt  
    … …  
    408406                  ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 
    409407                  ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8 
    410                   zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , epsi13 ) ) 
     408                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) ) 
    411409                  zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) )  
    412410                  zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
    … …  
    414412                  zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
    415413                     &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  
    416                   zds         = zfracs * sss_m(zji,zjj) - s_i_new(ji) 
    417                   s_i_new(ji) = zfracs * sss_m(zji,zjj) 
     414                  zds         = zfracs * sss_m(ii,ij) - s_i_new(ji) 
     415                  s_i_new(ji) = zfracs * sss_m(ii,ij) 
    418416               ENDIF ! fc_bo_i 
    419417            END DO ! ji 
    … …  
    432430                  &                            - rcp * ( ztmelts - rtt )                                    ) 
    433431               ! Basal growth rate = - F*dt / q 
    434                dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
     432               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
    435433               ! Salinity update 
    436434               ! entrapment during bottom growth 
    … …  
    453451            s_i_new(ji)   =  s_i_b(ji,nlay_i) 
    454452            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) 
    455             zfbase(ji)    =  zqfont_bo(ji) / rdt_ice     ! heat conservation test 
     453            zfbase(ji)    =  zqfont_bo(ji) * r1_rdtice     ! heat conservation test 
    456454            zdq_i(ji)     =  0._wp 
    457455            dh_i_bott(ji) =  0._wp 
    … …  
    461459      DO jk = nlay_i, 1, -1 
    462460         DO ji = kideb, kiut 
    463             IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN 
    464                ztmelts            =   - tmut * s_i_b(ji,jk) + rtt  
    465                IF( t_i_b(ji,jk) >= ztmelts ) THEN 
    466                   zdeltah(ji,jk)  = - zh_i(ji) 
    467                   dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk) 
    468                   zinnermelt(ji)   = 1._wp 
    469                ELSE  ! normal ablation 
    470                   zdeltah(ji,jk)  = - zqfont_bo(ji) / q_i_b(ji,jk) 
    471                   zqfont_bo(ji)   = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
    472                   zdeltah(ji,jk)  = MAX(zdeltah(ji,jk), - zh_i(ji) ) 
    473                   dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk) 
    474                   zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 
    475                   ! contribution to salt flux 
    476                   zji             = MOD( npb(ji) - 1, jpi ) + 1 
    477                   zjj             = ( npb(ji) - 1 ) / jpi + 1 
    478                   zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji)   ) * a_i_b(ji)   & 
    479                      &                              * MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice  
     461            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  >=  0._wp  ) THEN 
     462               ztmelts = - tmut * s_i_b(ji,jk) + rtt  
     463               IF( t_i_b(ji,jk) >= ztmelts ) THEN   !!gm : a comment is needed 
     464                  zdeltah   (ji,jk) = - zh_i(ji) 
     465                  dh_i_bott (ji   ) = dh_i_bott(ji) + zdeltah(ji,jk) 
     466                  zinnermelt(ji   ) = 1._wp 
     467               ELSE                                  ! normal ablation 
     468                  zdeltah  (ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) 
     469                  zqfont_bo(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
     470                  zdeltah  (ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) ) 
     471                  dh_i_bott(ji   ) = dh_i_bott(ji) + zdeltah(ji,jk) 
     472                  zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    480473               ENDIF 
     474               ! contribution to salt flux 
     475               zsfx_melt(ji) = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice  
    481476            ENDIF 
    482477         END DO ! ji 
    … …  
    493488               ENDIF 
    494489               IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN 
    495                   WRITE(numout,*) ' ALERTE heat loss for basal melt : zji, zjj, jl :', zji, zjj, jl 
     490                  WRITE(numout,*) ' ALERTE heat loss for basal melt : ii, ij, jl :', ii, ij, jl 
    496491                  WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji) 
    497492                  WRITE(numout,*) ' zfbase    : ', zfbase(ji) 
    … …  
    502497                  WRITE(numout,*) ' qlbbq_1d  : ', qlbbq_1d(ji) 
    503498                  WRITE(numout,*) ' s_i_new   : ', s_i_new(ji) 
    504                   WRITE(numout,*) ' sss_m     : ', sss_m(zji,zjj) 
     499                  WRITE(numout,*) ' sss_m     : ', sss_m(ii,ij) 
    505500                  WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 
    506501                  WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) ) 
    … …  
    531526         !                     ! excessive energy is sent to lateral ablation 
    532527         fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 )   & 
    533             &                          * ( zdhbf - dh_i_bott(ji) ) / rdt_ice 
     528            &                          * ( zdhbf - dh_i_bott(ji) ) * r1_rdtice 
    534529         dh_i_bott(ji)  = zdhbf 
    535530         !                     !since ice volume is only used for outputs, we keep it global for all categories 
    … …  
    538533         zhgnew   (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 
    539534         !                     ! diagnostic ( bottom ice growth ) 
    540          zji = MOD( npb(ji) - 1, jpi ) + 1 
    541          zjj = ( npb(ji) - 1 ) / jpi + 1 
    542          diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice 
    543          diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) / rdt_ice 
    544          diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice 
     535         ii = MOD( npb(ji) - 1, jpi ) + 1 
     536         ij = ( npb(ji) - 1 ) / jpi + 1 
     537         diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
     538         diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 
     539         diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
    545540      END DO 
    546541 
    … …  
    548543      ! 5.2 More than available ice melts 
    549544      !----------------------------------- 
    550       ! then heat applied minus heat content at previous time step 
    551       ! should equal heat remaining  
     545      ! then heat applied minus heat content at previous time step should equal heat remaining  
    552546      ! 
    553547      DO ji = kideb, kiut 
    554548         ! Adapt the remaining energy if too much ice melts 
    555549         !-------------------------------------------------- 
    556          zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 
     550         zihgnew =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice 
    557551         ! 0 if no more ice 
    558552         zhgnew    (ji) =         zihgnew   * zhgnew(ji)      ! ice thickness is put to 0 
    … …  
    562556         ! If snow remains, energy is used to melt snow 
    563557         zhni =  ht_s_b(ji)      ! snow depth at previous time step 
    564          zihg =  MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! 0 if snow  
     558         zihg =  MAX(  zzero , SIGN ( zone , - ht_s_b(ji) )  )   ! =0 if snow  
    565559 
    566560         ! energy of melting of remaining snow 
    567561         zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 ) 
    568562         zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 ) 
    569          zhnfi          =  zhni + zdhnm 
     563         zhnfi     =  zhni + zdhnm 
    570564         zfdt_final(ji) =  MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 ) 
    571565         ht_s_b(ji)     =  MAX( zzero , zhnfi ) 
    … …  
    581575         ! 
    582576         !                                              ! mass variation cumulated over category 
    583          rdmsnif_1d(ji) = rdmsnif_1d(ji) + zzfmass_s                     ! snow  
    584          rdmicif_1d(ji) = rdmicif_1d(ji) + zzfmass_i                     ! ice  
     577         rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow  
     578         rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice  
    585579 
    586580         ! Remaining heat to the ocean  
    587581         !--------------------------------- 
    588          focea(ji)  = - zfdt_final(ji) / rdt_ice         ! focea is in W.m-2 * dt 
    589  
    590       END DO 
    591  
    592       ftotal_fin (:) = zfdt_final(:)  / rdt_ice 
     582         focea(ji)  = - zfdt_final(ji) * r1_rdtice         ! focea is in W.m-2 * dt 
     583 
     584      END DO 
     585 
     586      ftotal_fin (:) = zfdt_final(:)  * r1_rdtice 
    593587 
    594588      !--------------------------- 
    … …  
    596590      !--------------------------- 
    597591      DO ji = kideb, kiut 
    598          zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   !1 if ice 
    599  
     592         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice 
     593         ! 
    600594         ! Salt flux 
    601          zji = MOD( npb(ji) - 1, jpi ) + 1 
    602          zjj = ( npb(ji) - 1 ) / jpi + 1 
    603          ! new lines 
    604          IF( num_sal == 4 ) THEN 
    605             fseqv_1d(ji) = fseqv_1d(ji) +        zihgnew  * zfsalt_melt(ji)                                & 
    606                &                        + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - bulk_sal   ) / rdt_ice 
    607          ELSE 
    608             fseqv_1d(ji) = fseqv_1d(ji) +        zihgnew  * zfsalt_melt(ji)                                & 
    609                &                        + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice 
    610          ENDIF 
     595         sfx_thd_1d(ji) = sfx_thd_1d(ji) +        zihgnew  * zsfx_melt(ji)               & 
     596            &                            - (1.0 - zihgnew) * zfmass_i (ji) * sm_i_b(ji)  * r1_rdtice 
     597         ! 
    611598         ! Heat flux 
    612599         ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1 
    613          ! excessive total ablation energy (focea) sent to the ocean 
     600         ! excessive total  ablation energy (focea) sent to the ocean 
    614601         qfvbq_1d(ji)  = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice 
    615602 
    616          zihic   = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) 
    617          ! equals 0 if ht_i = 0, 1 if ht_i gt 0 
     603         zihic   = 1.0 - MAX(  zzero , SIGN( zone , -ht_i_b(ji) )  )      ! equals 0 if ht_i = 0, 1 if ht_i gt 0 
    618604         fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji) 
    619          qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji)    * a_i_b(ji) * rdt_ice   & 
     605         qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea   (ji) * a_i_b(ji) * rdt_ice   & 
    620606            &                                    + ( 1.0 - zihic   ) * fscbq_1d(ji)             * rdt_ice 
    621607      END DO  ! ji 
    … …  
    656642         dmgwi_1d  (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 
    657643 
    658          rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic  
    659          rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn 
     644         ! All snow is thrown in the ocean, and seawater is taken to replace the volume 
     645         rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic * ( 1. - rhosn / rhoic ) 
     646         rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn 
    660647 
    661648         !        Equivalent salt flux (1) Snow-ice formation component 
    662649         !        ----------------------------------------------------- 
    663          zji = MOD( npb(ji) - 1, jpi ) + 1 
    664          zjj =    ( npb(ji) - 1 ) / jpi + 1 
    665  
    666          IF( num_sal /= 2 ) THEN   ;   zsm_snowice = sm_i_b(ji) 
    667          ELSE                      ;   zsm_snowice = ( rhoic - rhosn ) / rhoic * sss_m(zji,zjj)  
     650         ii = MOD( npb(ji) - 1, jpi ) + 1 
     651         ij =    ( npb(ji) - 1 ) / jpi + 1 
     652 
     653         IF( num_sal == 2 ) THEN   ;   zsm_snowice = sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic 
     654         ELSE                      ;   zsm_snowice = sm_i_b(ji)   
    668655         ENDIF 
    669          IF( num_sal == 4 ) THEN 
    670             fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal    ) * a_i_b(ji)   & 
    671                &                        * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 
    672          ELSE 
    673             fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - zsm_snowice ) * a_i_b(ji)   & 
    674                &                        * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 
    675          ENDIF 
     656         sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
     657         ! 
    676658         ! entrapment during snow ice formation 
    677659         i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 
    678660         isnowic      = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji)      ) ) * i_ice_switch 
    679          IF(  num_sal == 2  .OR.  num_sal == 4  )   & 
    680             dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) & 
    681             &               + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13)   & 
    682             &               - sm_i_b(ji) ) * isnowic      
     661         IF(  num_sal == 2  )   & 
     662            dsm_i_si_1d(ji) = (  zsm_snowice * dh_snowice(ji)    & 
     663            &                  + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 )   & 
     664            &                  - sm_i_b(ji) ) * isnowic      
    683665 
    684666         !  Actualize new snow and ice thickness. 
    … …  
    690672 
    691673         ! diagnostic ( snow ice growth ) 
    692          zji = MOD( npb(ji) - 1, jpi ) + 1 
    693          zjj =    ( npb(ji) - 1 ) / jpi + 1 
    694          diag_sni_gr(zji,zjj)  = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / rdt_ice 
     674         ii = MOD( npb(ji) - 1, jpi ) + 1 
     675         ij =    ( npb(ji) - 1 ) / jpi + 1 
     676         diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 
    695677         ! 
    696678      END DO !ji 
    697679      ! 
    698680      CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 
    699       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfsalt_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
     681      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zsfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
    700682      CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) 
    701683      CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r3610 r3625  
    1515   !!   'key_lim3'                                      LIM3 sea-ice model 
    1616   !!---------------------------------------------------------------------- 
    17    USE par_oce          ! ocean parameters 
    18    USE phycst           ! physical constants (ocean directory)  
    19    USE ice              ! LIM-3 variables 
    20    USE par_ice          ! LIM-3 parameters 
    21    USE thd_ice          ! LIM-3: thermodynamics 
    22    USE in_out_manager   ! I/O manager 
    23    USE lib_mpp          ! MPP library 
    24    USE wrk_nemo         ! work arrays 
     17   USE par_oce        ! ocean parameters 
     18   USE phycst         ! physical constants (ocean directory)  
     19   USE ice            ! LIM-3 variables 
     20   USE par_ice        ! LIM-3 parameters 
     21   USE thd_ice        ! LIM-3: thermodynamics 
     22   USE in_out_manager ! I/O manager 
     23   USE lib_mpp        ! MPP library 
     24   USE wrk_nemo       ! work arrays 
     25   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2526 
    2627   IMPLICIT NONE 
    … …  
    3334 
    3435   !!---------------------------------------------------------------------- 
    35    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     36   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    3637   !! $Id$ 
    3738   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    147148      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis 
    148149      REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms 
    149       !!------------------------------------------------------------------ 
    150        
     150      !!------------------------------------------------------------------      
    151151      !  
    152152      !------------------------------------------------------------------------------! 
    … …  
    156156      DO ji = kideb , kiut 
    157157         ! is there snow or not 
    158          isnow(ji)= INT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  ) 
     158         isnow(ji)= INT(  1._wp - MAX(  0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  ) 
    159159         ! surface temperature of fusion 
    160160!!gm ???  ztfs(ji) = rtt !!!???? 
    … …  
    201201      DO ji = kideb , kiut 
    202202         ! switches 
    203          isnow(ji) = INT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  )  
     203         isnow(ji) = INT(  1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  )  
    204204         ! hs > 0, isnow = 1 
    205205         zhsu (ji) = hnzst  ! threshold for the computation of i0 
    … …  
    262262      ! just to check energy conservation 
    263263      DO ji = kideb, kiut 
    264          ii                = MOD( npb(ji) - 1, jpi ) + 1 
    265          ij                = ( npb(ji) - 1 ) / jpi + 1 
     264         ii = MOD( npb(ji) - 1 , jpi ) + 1 
     265         ij =    ( npb(ji) - 1 ) / jpi + 1 
    266266         fstroc(ii,ij,jl) = zradtr_i(ji,nlay_i) 
    267267      END DO 
    … …  
    273273         END DO 
    274274      END DO 
    275  
    276275 
    277276      ! 
    … …  
    662661 
    663662            ! surface temperature 
    664             isnow(ji)     = INT(1.0-max(0.0,sign(1.0,-ht_s_b(ji)))) 
     663            isnow(ji)     = INT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) )  )  ) 
    665664            ztsuoldit(ji) = t_su_b(ji) 
    666             IF (t_su_b(ji) .LT. ztfs(ji)) & 
     665            IF( t_su_b(ji) < ztfs(ji) )  & 
    667666               t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( isnow(ji)*t_s_b(ji,1)   & 
    668667               &          + (1.0-isnow(ji))*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r3294 r3625  
    1616   !!   'key_lim3'                                      LIM3 sea-ice model 
    1717   !!---------------------------------------------------------------------- 
    18    !!   lim_thd_ent : ice redistribution of enthalpy 
     18   !!   lim_thd_ent   : ice redistribution of enthalpy 
    1919   !!---------------------------------------------------------------------- 
    20    USE par_oce          ! ocean parameters 
    21    USE dom_oce          ! domain variables 
    22    USE domain           ! 
    23    USE phycst           ! physical constants 
    24    USE ice              ! LIM variables 
    25    USE par_ice          ! LIM parameters 
    26    USE thd_ice          ! LIM thermodynamics 
    27    USE limvar           ! LIM variables 
    28    USE in_out_manager   ! I/O manager 
    29    USE lib_mpp          ! MPP library 
    30    USE wrk_nemo         ! work arrays 
     20   USE par_oce        ! ocean parameters 
     21   USE dom_oce        ! domain variables 
     22   USE domain         ! 
     23   USE phycst         ! physical constants 
     24   USE ice            ! LIM variables 
     25   USE par_ice        ! LIM parameters 
     26   USE thd_ice        ! LIM thermodynamics 
     27   USE limvar         ! LIM variables 
     28   USE in_out_manager ! I/O manager 
     29   USE lib_mpp        ! MPP library 
     30   USE wrk_nemo       ! work arrays 
     31   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3132 
    3233   IMPLICIT NONE 
    … …  
    4344 
    4445   !!---------------------------------------------------------------------- 
    45    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     46   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    4647   !! $Id$ 
    4748   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    408409      IF ( con_i ) THEN 
    409410         DO ji = kideb, kiut 
    410             IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN 
     411            IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice  > 1.0e-6 ) THEN 
    411412               zji                 = MOD( npb(ji) - 1, jpi ) + 1 
    412413               zjj                 = ( npb(ji) - 1 ) / jpi + 1 
    413414               WRITE(numout,*) ' violation of heat conservation : ',             & 
    414                   ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice 
     415                  ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 
    415416               WRITE(numout,*) ' ji, jj   : ', zji, zjj 
    416417               WRITE(numout,*) ' ht_s_b   : ', ht_s_b(ji) 
    417                WRITE(numout,*) ' zqts_in  : ', zqts_in(ji) / rdt_ice 
    418                WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) / rdt_ice 
     418               WRITE(numout,*) ' zqts_in  : ', zqts_in (ji) * r1_rdtice 
     419               WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) * r1_rdtice 
    419420               WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji) 
    420421               WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji) 
    … …  
    526527         ! bottom formation temperature 
    527528         ztform = t_i_b(ji,nlay_i) 
    528          IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) ztform = t_bo_b(ji) 
     529         IF(  num_sal == 2  )  ztform = t_bo_b(ji) 
    529530         qm0(ji,nbot0(ji)) = ( 1.0 - icboswi(ji) )*qm0(ji,nbot0(ji))             &   ! case of melting ice 
    530531            &              + icboswi(ji) * rhoic * ( cpic*(ztmelts-ztform)       &   ! case of forming ice 
    … …  
    622623      ! 
    623624      DO ji = kideb, kiut 
    624          IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN 
     625         IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice  > 1.0e-6 ) THEN 
    625626            zji                 = MOD( npb(ji) - 1, jpi ) + 1 
    626627            zjj                 = ( npb(ji) - 1 ) / jpi + 1 
    627             WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice 
     628            WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 
    628629            WRITE(numout,*) ' ji, jj   : ', zji, zjj 
    629630            WRITE(numout,*) ' ht_i_b   : ', ht_i_b(ji) 
    630             WRITE(numout,*) ' zqti_in  : ', zqti_in(ji) / rdt_ice 
    631             WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) / rdt_ice 
     631            WRITE(numout,*) ' zqti_in  : ', zqti_in (ji) * r1_rdtice 
     632            WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice 
    632633            WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) 
    633634            WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r3294 r3625  
    1313   !!   'key_lim3'                                      LIM3 sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    15    !!   lim_lat_acr    : lateral accretion of ice 
    16    !!---------------------------------------------------------------------- 
    17    USE par_oce          ! ocean parameters 
    18    USE dom_oce          ! domain variables 
    19    USE phycst           ! physical constants 
    20    USE sbc_oce          ! Surface boundary condition: ocean fields 
    21    USE sbc_ice          ! Surface boundary condition: ice fields 
    22    USE thd_ice          ! LIM thermodynamics 
    23    USE dom_ice          ! LIM domain 
    24    USE par_ice          ! LIM parameters 
    25    USE ice              ! LIM variables 
    26    USE limtab           ! LIM 2D <==> 1D 
    27    USE limcons          ! LIM conservation 
    28    USE in_out_manager   ! I/O manager 
    29    USE lib_mpp          ! MPP library 
    30    USE wrk_nemo         ! work arrays 
     15   !!   lim_lat_acr   : lateral accretion of ice 
     16   !!---------------------------------------------------------------------- 
     17   USE par_oce        ! ocean parameters 
     18   USE dom_oce        ! domain variables 
     19   USE phycst         ! physical constants 
     20   USE sbc_oce        ! Surface boundary condition: ocean fields 
     21   USE sbc_ice        ! Surface boundary condition: ice fields 
     22   USE thd_ice        ! LIM thermodynamics 
     23   USE dom_ice        ! LIM domain 
     24   USE par_ice        ! LIM parameters 
     25   USE ice            ! LIM variables 
     26   USE limtab         ! LIM 2D <==> 1D 
     27   USE limcons        ! LIM conservation 
     28   USE in_out_manager ! I/O manager 
     29   USE lib_mpp        ! MPP library 
     30   USE wrk_nemo       ! work arrays 
     31   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3132 
    3233   IMPLICIT NONE 
    … …  
    4546 
    4647   !!---------------------------------------------------------------------- 
    47    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     48   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    4849   !! $Id$ 
    4950   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    7778      !!               update ht_s_b, ht_i_b and tbif_1d(:,:)       
    7879      !!------------------------------------------------------------------------ 
    79       INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices 
    80       INTEGER ::   layer, nbpac     ! local integers  
    81       INTEGER ::   zji, zjj, iter   !   -       - 
    82       REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde  ! local scalars 
     80      INTEGER  ::   ji,jj,jk,jl,jm   ! dummy loop indices 
     81      INTEGER  ::   layer, nbpac     ! local integers  
     82      INTEGER  ::   zji, zjj, iter   !   -       - 
     83      REAL(wp) ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde   ! local scalars 
    8384      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      - 
    8485      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      - 
     86      REAL(wp) ::   zcoef                                                       !   -      - 
    8587      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness 
    8688      CHARACTER (len = 15) :: fieldid 
    … …  
    143145      ! 1) Conservation check and changes in each ice category 
    144146      !------------------------------------------------------------------------------! 
    145       IF ( con_i ) THEN 
    146          CALL lim_column_sum (jpl, v_i, vt_i_init) 
    147          CALL lim_column_sum (jpl, v_s, vt_s_init) 
    148          CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init) 
    149          CALL lim_column_sum (jpl,  e_s(:,:,1,:) , et_s_init) 
     147      IF( con_i ) THEN 
     148         CALL lim_column_sum        ( jpl, v_i          , vt_i_init) 
     149         CALL lim_column_sum        ( jpl, v_s          , vt_s_init) 
     150         CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init) 
     151         CALL lim_column_sum        ( jpl, e_s(:,:,1,:) , et_s_init) 
    150152      ENDIF 
    151153 
    … …  
    158160               DO ji = 1, jpi 
    159161                  !Energy of melting q(S,T) [J.m-3] 
    160                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / & 
    161                      MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * & 
    162                      nlay_i 
    163                   zindb      = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl))) !0 if no ice and 1 if yes 
    164                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl)*unit_fac*zindb 
     162                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * nlay_i 
     163                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) )  )   !0 if no ice and 1 if yes 
     164                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 
    165165               END DO 
    166166            END DO 
    … …  
    182182      !  
    183183 
    184       zvrel(:,:) = 0.0 
     184      zvrel(:,:) = 0._wp 
    185185 
    186186      ! Default new ice thickness  
    187       DO jj = 1, jpj 
    188          DO ji = 1, jpi 
    189             hicol(ji,jj) = hiccrit(1) 
    190          END DO 
    191       END DO 
    192  
    193       IF (fraz_swi.eq.1.0) THEN 
     187      hicol(:,:) = hiccrit(1) 
     188 
     189      IF( fraz_swi == 1._wp ) THEN 
    194190 
    195191         !-------------------- 
    196192         ! Physical constants 
    197193         !-------------------- 
    198          hicol(:,:) = 0.0 
     194         hicol(:,:) = 0._wp 
    199195 
    200196         zhicrit = 0.04 ! frazil ice thickness 
    … …  
    211207                  !------------- 
    212208                  ! C-grid wind stress components 
    213                   ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  ) & 
    214                      &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) / 2.0 
    215                   ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1) & 
    216                      &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) / 2.0 
     209                  ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  )   & 
     210                     &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) * 0.5_wp 
     211                  ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1)   & 
     212                     &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) * 0.5_wp 
    217213                  ! Square root of wind stress 
    218214                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
    … …  
    228224                  !------------------- 
    229225                  ! C-grid ice velocity 
    230                   zindb = MAX(0.0, SIGN(1.0, at_i(ji,jj) )) 
    231                   zvgx  = zindb * ( u_ice(ji-1,jj  ) * tmu(ji-1,jj  ) & 
    232                      + u_ice(ji,jj    ) * tmu(ji  ,jj  ) ) / 2.0 
    233                   zvgy  = zindb * ( v_ice(ji  ,jj-1) * tmv(ji  ,jj-1) & 
    234                      + v_ice(ji,jj    ) * tmv(ji  ,jj  ) ) / 2.0 
     226                  zindb = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  ) 
     227                  zvgx  = zindb * (  u_ice(ji-1,jj  ) * tmu(ji-1,jj  )    & 
     228                     &             + u_ice(ji,jj    ) * tmu(ji  ,jj  )  ) * 0.5_wp 
     229                  zvgy  = zindb * (  v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)    & 
     230                     &             + v_ice(ji,jj    ) * tmv(ji  ,jj  )  ) * 0.5_wp 
    235231 
    236232                  !----------------------------------- 
    … …  
    238234                  !----------------------------------- 
    239235                  ! absolute relative velocity 
    240                   zvrel2        = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) + & 
    241                      ( zvfry - zvgy ) * ( zvfry - zvgy )   & 
    242                      , 0.15 * 0.15 ) 
    243                   zvrel(ji,jj)  = SQRT(zvrel2) 
     236                  zvrel2 = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   & 
     237                     &         + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) 
     238                  zvrel(ji,jj)  = SQRT( zvrel2 ) 
    244239 
    245240                  !--------------------- 
    … …  
    247242                  !--------------------- 
    248243                  hicol(ji,jj) = zhicrit + 0.1  
    249                   hicol(ji,jj) = zhicrit + hicol(ji,jj) /      &  
    250                      ( hicol(ji,jj) * hicol(ji,jj) - & 
    251                      zhicrit * zhicrit ) * ztwogp * zvrel2 
     244                  hicol(ji,jj) = zhicrit +   hicol(ji,jj)    & 
     245                     &                   / ( hicol(ji,jj) * hicol(ji,jj) -  zhicrit * zhicrit ) * ztwogp * zvrel2 
     246 
     247!!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1) 
     248!!gm                                                   = zhicrit**2 + 0.2*zhicrit +0.01 
     249!!gm                therefore the 2 lines with hicol can be replaced by 1 line: 
     250!!gm              hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2 
     251!!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop 
    252252 
    253253                  iter = 1 
    … …  
    284284      DO jj = 1, jpj 
    285285         DO ji = 1, jpi 
    286             IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN 
     286            IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN 
    287287               nbpac = nbpac + 1 
    288288               npac( nbpac ) = (jj - 1) * jpi + ji 
    289                IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN 
    290                   jiindex_1d = nbpac 
    291                ENDIF 
     289               IF( ji == jiindx .AND. jj == jjindx )   jiindex_1d = nbpac 
    292290            ENDIF 
    293291         END DO 
    294292      END DO 
    295293 
    296       IF( ln_nicep ) THEN 
    297          WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 
    298       ENDIF 
     294      IF( ln_nicep )   WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 
    299295 
    300296      !------------------------------ 
    … …  
    306302      IF ( nbpac > 0 ) THEN 
    307303 
    308          CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         ,       & 
    309             jpi, jpj, npac(1:nbpac) ) 
     304         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) ) 
    310305         DO jl = 1, jpl 
    311             CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl)  , a_i(:,:,jl)  ,       & 
    312                jpi, jpj, npac(1:nbpac) ) 
    313             CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl)  , v_i(:,:,jl)  ,       & 
    314                jpi, jpj, npac(1:nbpac) ) 
    315             CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) ,       & 
    316                jpi, jpj, npac(1:nbpac) ) 
    317             CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl),       & 
    318                jpi, jpj, npac(1:nbpac) ) 
     306            CALL tab_2d_1d( nbpac, za_i_ac  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
     307            CALL tab_2d_1d( nbpac, zv_i_ac  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
     308            CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
     309            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    319310            DO jk = 1, nlay_i 
    320                CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , & 
    321                   jpi, jpj, npac(1:nbpac) ) 
     311               CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 
    322312            END DO ! jk 
    323313         END DO ! jl 
    324314 
    325          CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif ,              & 
    326             jpi, jpj, npac(1:nbpac) ) 
    327          CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif ,              & 
    328             jpi, jpj, npac(1:nbpac) ) 
    329          CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo  ,              & 
    330             jpi, jpj, npac(1:nbpac) ) 
    331          CALL tab_2d_1d( nbpac, fseqv_1d  (1:nbpac)     , fseqv ,              & 
    332             jpi, jpj, npac(1:nbpac) ) 
    333          CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol ,              & 
    334             jpi, jpj, npac(1:nbpac) ) 
    335          CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel ,              & 
    336             jpi, jpj, npac(1:nbpac) ) 
     315         CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif  , jpi, jpj, npac(1:nbpac) ) 
     316         CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif  , jpi, jpj, npac(1:nbpac) ) 
     317         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) ) 
     318         CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac)     , sfx_thd, jpi, jpj, npac(1:nbpac) ) 
     319         CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice, jpi, jpj, npac(1:nbpac) ) 
     320         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) ) 
     321         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) ) 
    337322 
    338323         !------------------------------------------------------------------------------! 
    … …  
    344329         !---------------------- 
    345330         DO ji = 1, nbpac 
    346             zh_newice(ji)     = hiccrit(1) 
    347          END DO 
    348          IF ( fraz_swi .EQ. 1.0 ) zh_newice(:) = hicol_b(:) 
     331            zh_newice(ji) = hiccrit(1) 
     332         END DO 
     333         IF( fraz_swi == 1.0 )  zh_newice(:) = hicol_b(:) 
    349334 
    350335         !---------------------- 
    … …  
    352337         !---------------------- 
    353338 
    354          IF ( num_sal .EQ. 1 ) THEN 
    355             zs_newice(:)      =   bulk_sal 
    356          ENDIF ! num_sal 
    357  
    358          IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN 
    359  
    360             DO ji = 1, nbpac 
    361                zs_newice(ji)  =   MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max ) 
    362                zji            =   MOD( npac(ji) - 1, jpi ) + 1 
    363                zjj            =   ( npac(ji) - 1 ) / jpi + 1 
    364                zs_newice(ji)  =   MIN( 0.5*sss_m(zji,zjj) , zs_newice(ji) ) 
    365             END DO ! jl 
    366  
    367          ENDIF ! num_sal 
    368  
    369          IF ( num_sal .EQ. 3 ) THEN 
    370             zs_newice(:)      =   2.3 
    371          ENDIF ! num_sal 
     339         SELECT CASE ( num_sal ) 
     340         CASE ( 1 )                    ! Sice = constant  
     341            zs_newice(:) = bulk_sal 
     342         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)] 
     343            DO ji = 1, nbpac 
     344               zji =   MOD( npac(ji) - 1 , jpi ) + 1 
     345               zjj =      ( npac(ji) - 1 ) / jpi + 1 
     346               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(zji,zjj)  ) 
     347            END DO 
     348         CASE ( 3 )                    ! Sice = F(z) [multiyear ice] 
     349            zs_newice(:) =   2.3 
     350         END SELECT 
     351 
    372352 
    373353         !------------------------- 
    … …  
    376356         ! We assume that new ice is formed at the seawater freezing point 
    377357         DO ji = 1, nbpac 
    378             ztmelts           = - tmut * zs_newice(ji) + rtt ! Melting point (K) 
    379             ze_newice(ji)     =   rhoic * ( cpic * ( ztmelts - t_bo_b(ji) )    & 
    380                + lfus * ( 1.0 - ( ztmelts - rtt )   & 
    381                / ( t_bo_b(ji) - rtt ) )           & 
    382                - rcp * ( ztmelts-rtt ) ) 
    383             ze_newice(ji)     =   MAX( ze_newice(ji) , 0.0 ) +                 & 
    384                MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) )   &  
    385                * rhoic * lfus 
     358            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K) 
     359            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             & 
     360               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   & 
     361               &                       - rcp  *         ( ztmelts - rtt )  ) 
     362            ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    & 
     363               &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus 
    386364         END DO ! ji 
    387365         !---------------- 
    … …  
    389367         !---------------- 
    390368         DO ji = 1, nbpac 
    391             zo_newice(ji)     = 0.0 
     369            zo_newice(ji) = 0._wp 
    392370         END DO ! ji 
    393371 
    … …  
    396374         !-------------------------- 
    397375         DO ji = 1, nbpac 
    398             zqbgow(ji)        = qldif_1d(ji) - qcmif_1d(ji) !<0 
     376            zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)    !<0 
    399377         END DO ! ji 
    400378 
    … …  
    403381         !------------------- 
    404382         DO ji = 1, nbpac 
    405             zv_newice(ji)     = - zqbgow(ji) / ze_newice(ji) 
     383            zv_newice(ji) = - zqbgow(ji) / ze_newice(ji) 
    406384 
    407385            ! A fraction zfrazb of frazil ice is accreted at the ice bottom 
    408             zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) )     &  
    409                + 1.0 ) / 2.0 * maxfrazb 
    410             zdh_frazb(ji) = zfrazb*zv_newice(ji) 
     386            zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 
     387            zdh_frazb(ji) =         zfrazb   * zv_newice(ji) 
    411388            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 
    412389         END DO 
    … …  
    415392         ! Salt flux due to new ice growth 
    416393         !--------------------------------- 
    417          IF ( ( num_sal .EQ. 4 ) ) THEN  
    418             DO ji = 1, nbpac 
    419                zji            = MOD( npac(ji) - 1, jpi ) + 1 
    420                zjj            = ( npac(ji) - 1 ) / jpi + 1 
    421                fseqv_1d(ji)   = fseqv_1d(ji) +                                     & 
    422                   ( sss_m(zji,zjj) - bulk_sal      ) * rhoic *       & 
    423                   zv_newice(ji) / rdt_ice 
    424             END DO 
    425          ELSE 
    426             DO ji = 1, nbpac 
    427                zji            = MOD( npac(ji) - 1, jpi ) + 1 
    428                zjj            = ( npac(ji) - 1 ) / jpi + 1 
    429                fseqv_1d(ji)   = fseqv_1d(ji) +                                     & 
    430                   ( sss_m(zji,zjj) - zs_newice(ji) ) * rhoic *       & 
    431                   zv_newice(ji) / rdt_ice 
    432             END DO ! ji 
    433          ENDIF 
     394         ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine) 
     395         DO ji = 1, nbpac 
     396            sfx_thd_1d(ji) = sfx_thd_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice 
     397            rdm_ice_1d(ji) = rdm_ice_1d(ji) +                 rhoic * zv_newice(ji) 
     398         END DO ! ji 
    434399 
    435400         !------------------------------------ 
    … …  
    437402         !------------------------------------ 
    438403         DO ji = 1, nbpac 
    439             ! Volume 
    440             zji                  = MOD( npac(ji) - 1, jpi ) + 1 
    441             zjj                  = ( npac(ji) - 1 ) / jpi + 1 
    442             vt_i_init(zji,zjj)   = vt_i_init(zji,zjj) + zv_newice(ji) 
    443             ! Energy 
    444             zde                  = ze_newice(ji) / unit_fac 
    445             zde                  = zde * area(zji,zjj) * zv_newice(ji) 
    446             et_i_init(zji,zjj)   = et_i_init(zji,zjj) + zde 
     404            zji = MOD( npac(ji) - 1 , jpi ) + 1 
     405            zjj =    ( npac(ji) - 1 ) / jpi + 1 
     406            ! 
     407            zde = ze_newice(ji) / unit_fac * area(zji,zjj) * zv_newice(ji) 
     408            ! 
     409            vt_i_init(zji,zjj) = vt_i_init(zji,zjj) + zv_newice(ji)             ! volume 
     410            et_i_init(zji,zjj) = et_i_init(zji,zjj) + zde                       ! Energy 
     411 
    447412         END DO 
    448413 
    449414         ! keep new ice volume in memory 
    450          CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , & 
    451             jpi, jpj ) 
     415         CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj ) 
    452416 
    453417         !----------------- 
    … …  
    455419         !----------------- 
    456420         DO ji = 1, nbpac 
    457             za_newice(ji)     = zv_newice(ji) / zh_newice(ji) 
    458             ! diagnostic 
    459             zji                  = MOD( npac(ji) - 1, jpi ) + 1 
    460             zjj                  = ( npac(ji) - 1 ) / jpi + 1 
    461     &nb