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 for branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Ignore:
Timestamp:
2012-11-21T14:19:18+01:00 (11 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/NEMOGCM/NEMO/LIM_SRC_3
Files:
28 edited

Legend:

Unmodified
Added
Removed

  • 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             diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice 
     421            zji = MOD( npac(ji) - 1 , jpi ) + 1 
     422            zjj =    ( npac(ji) - 1 ) / jpi + 1 
     423            za_newice(ji) = zv_newice(ji) / zh_newice(ji) 
     424            diag_lat_gr(zji,zjj) = zv_newice(ji) * r1_rdtice 
    462425         END DO !ji 
    463426 
    … …  
    476439         !------------------------------------------- 
    477440         ! If lateral ice growth gives an ice concentration gt 1, then 
    478          ! we keep the excessive volume in memory and attribute it later 
    479          ! to bottom accretion 
    480          DO ji = 1, nbpac 
    481             ! vectorize 
    482             IF ( za_newice(ji) .GT. ( 1.0 - zat_i_ac(ji) ) ) THEN 
    483                zda_res(ji)    = za_newice(ji) - (1.0 - zat_i_ac(ji) ) 
    484                zdv_res(ji)    = zda_res(ji) * zh_newice(ji)  
    485                za_newice(ji)  = za_newice(ji) - zda_res(ji) 
    486                zv_newice(ji)  = zv_newice(ji) - zdv_res(ji) 
     441         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
     442         DO ji = 1, nbpac 
     443            IF ( za_newice(ji)  >  ( 1._wp - zat_i_ac(ji) ) ) THEN 
     444               zda_res(ji)   = za_newice(ji) - (1.0 - zat_i_ac(ji) ) 
     445               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji)  
     446               za_newice(ji) = za_newice(ji) - zda_res  (ji) 
     447               zv_newice(ji) = zv_newice(ji) - zdv_res  (ji) 
    487448            ELSE 
    488                zda_res(ji) = 0.0 
    489                zdv_res(ji) = 0.0 
     449               zda_res(ji) = 0._wp 
     450               zdv_res(ji) = 0._wp 
    490451            ENDIF 
    491452         END DO ! ji 
    … …  
    497458         DO jl = 1, jpl 
    498459            DO ji = 1, nbpac 
    499                IF(  hi_max   (jl-1)  <  zh_newice(ji)   .AND.   & 
    500                   & zh_newice(ji)    <= hi_max   (jl)         ) THEN 
     460               IF(  hi_max   (jl-1)  <   zh_newice(ji)   .AND.   & 
     461                  & zh_newice(ji)    <=  hi_max   (jl)         ) THEN 
    501462                  za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji) 
    502463                  zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji) 
    … …  
    504465                  zcatac  (ji)    = jl 
    505466               ENDIF 
    506             END DO ! ji 
    507          END DO ! jl 
     467            END DO 
     468         END DO 
    508469 
    509470         !---------------------------------- 
    … …  
    521482            DO ji = 1, nbpac 
    522483               jl = zcatac(ji) 
    523                zqold   = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 
     484               zqold   = ze_i_ac(ji,jk,jl)      ! [ J.m-3 ] 
    524485               zalphai = MIN( zhice_old(ji,jl) *   jk       / nlay_i , zh_newice(ji) )   & 
    525486                  &    - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i , zh_newice(ji) ) 
    … …  
    527488                  + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / nlay_i   & 
    528489                  + za_newice(ji)  * ze_newice(ji) * zalphai                                       & 
    529                   + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( ( zv_i_ac(ji,jl) ) / nlay_i ) 
     490                  + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( zv_i_ac(ji,jl) / nlay_i ) 
    530491            END DO 
    531492         END DO 
    … …  
    567528               zdhicbot (ji,jl) = zdv_res(ji)    / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    & 
    568529                  &             +  zindb * zdh_frazb(ji)                               ! frazil ice may coalesce 
    569                zdummy(ji,jl)    = zv_i_ac(ji,jl)/MAX(za_i_ac(ji,jl),epsi10)*zindb      ! thickness of residual ice 
     530               zdummy(ji,jl)    = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb      ! thickness of residual ice 
    570531            END DO 
    571532         END DO 
    … …  
    628589         ! Update salinity 
    629590         !----------------- 
    630          IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN 
     591         IF(  num_sal == 2  ) THEN      ! Sice = F(z,t) 
    631592            DO jl = 1, jpl 
    632593               DO ji = 1, nbpac 
    633                   zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )  ! 0 if no ice and 1 if yes 
     594                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )   ! 0 if no ice and 1 if yes 
    634595                  zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl) 
    635596                  zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb 
    … …  
    645606            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 
    646607            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 
    647             IF (  num_sal == 2  .OR.  num_sal == 4  )   & 
     608            IF (  num_sal == 2  )   & 
    648609               CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 
    649610            DO jk = 1, nlay_i 
    … …  
    651612            END DO 
    652613         END DO 
    653          CALL tab_1d_2d( nbpac, fseqv , npac(1:nbpac), fseqv_1d  (1:nbpac) , jpi, jpj ) 
     614         CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj ) 
     615         CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj ) 
    654616         ! 
    655617      ENDIF ! nbpac > 0 
    … …  
    660622      DO jl = 1, jpl 
    661623         DO jk = 1, nlay_i          ! heat content in 10^9 Joules 
    662             e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i  / unit_fac  
     624            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i / unit_fac  
    663625         END DO 
    664626      END DO 

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

    r3294 r3625  
    1212   !!   'key_lim3'                                      LIM-3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    !!   lim_thd_sal : salinity variations in the ice 
    15    !!---------------------------------------------------------------------- 
    16    USE par_oce          ! ocean parameters 
    17    USE phycst           ! physical constants (ocean directory) 
    18    USE sbc_oce          ! Surface boundary condition: ocean fields 
    19    USE ice              ! LIM variables 
    20    USE par_ice          ! LIM parameters 
    21    USE thd_ice          ! LIM thermodynamics 
    22    USE limvar           ! LIM variables 
    23    USE in_out_manager   ! I/O manager 
    24    USE lib_mpp          ! MPP library 
    25    USE wrk_nemo         ! work arrays 
     14   !!   lim_thd_sal   : salinity variations in the ice 
     15   !!---------------------------------------------------------------------- 
     16   USE par_oce        ! ocean parameters 
     17   USE phycst         ! physical constants (ocean directory) 
     18   USE sbc_oce        ! Surface boundary condition: ocean fields 
     19   USE ice            ! LIM variables 
     20   USE par_ice        ! LIM parameters 
     21   USE thd_ice        ! LIM thermodynamics 
     22   USE limvar         ! LIM variables 
     23   USE in_out_manager ! I/O manager 
     24   USE lib_mpp        ! MPP library 
     25   USE wrk_nemo       ! work arrays 
     26   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2627 
    2728   IMPLICIT NONE 
    … …  
    3233 
    3334   !!---------------------------------------------------------------------- 
    34    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     35   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    3536   !! $Id$ 
    3637   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    4445      !! ** Purpose :   computes new salinities in the ice 
    4546      !! 
    46       !! ** Method  :  4 possibilities 
    47       !!               -> num_sal = 1 -> constant salinity for z,t 
    48       !!               -> num_sal = 2 -> S = S(z,t) [simple Vancoppenolle et al 2005] 
    49       !!               -> num_sal = 3 -> S = S(z)   [multiyear ice] 
    50       !!               -> num_sal = 4 -> S = S(h)   [Cox and Weeks 74] 
     47      !! ** Method  :  3 possibilities 
     48      !!               -> num_sal = 1 -> Sice = cst    [ice salinity constant in both time & space]  
     49      !!               -> num_sal = 2 -> Sice = S(z,t) [Vancoppenolle et al. 2005] 
     50      !!               -> num_sal = 3 -> Sice = S(z)   [multiyear ice] 
    5151      !!--------------------------------------------------------------------- 
    52       INTEGER, INTENT(in) ::  kideb, kiut   ! thickness category index 
     52      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index 
    5353      ! 
    5454      INTEGER  ::   ji, jk     ! dummy loop indices  
    55       INTEGER  ::   zji, zjj   ! local integers 
    5655      REAL(wp) ::   zsold, iflush, iaccrbo, igravdr, isnowic, i_ice_switch,  ztmelts   ! local scalars 
    5756      REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars 
    … …  
    6463      ! 1) Constant salinity, constant in time                                       | 
    6564      !------------------------------------------------------------------------------| 
    66 !!gm comment: if num_sal = 1 s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !! 
    67       IF( num_sal == 1 ) THEN 
    68          ! 
    69          DO jk = 1, nlay_i 
    70             DO ji = kideb, kiut 
    71                s_i_b(ji,jk) =  bulk_sal 
    72             END DO ! ji 
    73          END DO ! jk 
    74          ! 
    75          DO ji = kideb, kiut 
    76             sm_i_b(ji)      =  bulk_sal  
    77          END DO ! ji 
    78          ! 
     65!!gm comment: if num_sal = 1 s_i_new, s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !! 
     66!!gm           ===>>>   simplification of almost all test on num_sal value 
     67      IF(  num_sal == 1  ) THEN 
     68            s_i_b  (kideb:kiut,1:nlay_i) =  bulk_sal 
     69            sm_i_b (kideb:kiut)          =  bulk_sal  
     70            s_i_new(kideb:kiut)          =  bulk_sal 
    7971      ENDIF 
    8072 
    … …  
    8375      !------------------------------------------------------------------------------| 
    8476 
    85       IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN 
     77      IF(  num_sal == 2  ) THEN 
    8678 
    8779         !--------------------------------- 
    … …  
    118110            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice  
    119111            !                                   ! drainage by flushing   
    120             dsm_i_fl_1d(ji) = - iflush * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice 
     112            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice 
    121113 
    122114            !----------------- 
    … …  
    133125         END DO ! ji 
    134126 
    135          ! Salinity profile 
    136          CALL lim_var_salprof1d( kideb, kiut ) 
     127         CALL lim_var_salprof1d( kideb, kiut )         ! Salinity profile 
     128 
    137129 
    138130         !---------------------------- 
    … …  
    143135!!gm useless 
    144136            ! iflush  : 1 if summer  
    145             iflush  =  MAX( 0._wp , SIGN ( 1._wp , t_su_b(ji) - rtt ) )  
     137            iflush  =  MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt        ) )  
    146138            ! igravdr : 1 if t_su lt t_bo 
    147             igravdr =  MAX( 0._wp , SIGN ( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )  
     139            igravdr =  MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )  
    148140            ! iaccrbo : 1 if bottom accretion 
    149             iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) ) 
     141            iaccrbo =  MAX( 0._wp , SIGN( 1._wp , dh_i_bott(ji)           ) ) 
    150142!!gm end useless 
    151143            ! 
    … …  
    157149         !---------------------------- 
    158150         DO ji = kideb, kiut 
    159             i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) ) 
    160             fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji)         & 
    161                &         * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice 
    162             IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp 
    163          END DO ! ji 
     151            i_ice_switch = 1._wp - MAX(  0._wp, SIGN( 1._wp , - ht_i_b(ji) )  ) 
     152            sfx_bri_1d(ji) = sfx_bri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji)         & 
     153               &           * ( MAX( dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji) , sm_i_b(ji) - zsiold(ji) ) ) * r1_rdtice 
     154         END DO 
    164155 
    165156         ! Only necessary for conservation check since salinity is modified 
    … …  
    179170         END DO 
    180171         ! 
    181       ENDIF ! num_sal .EQ. 2 
     172      ENDIF  
    182173 
    183174      !------------------------------------------------------------------------------| 
    … …  
    185176      !------------------------------------------------------------------------------| 
    186177 
    187       IF( num_sal == 3 )   CALL lim_var_salprof1d( kideb, kiut ) 
    188  
    189       !------------------------------------------------------------------------------| 
    190       !  Module 4 : Constant salinity varying in time                                | 
    191       !------------------------------------------------------------------------------| 
    192  
    193       IF( num_sal == 5 ) THEN      ! Cox and Weeks, 1974 
    194          ! 
    195          DO ji = kideb, kiut 
    196             zsold = sm_i_b(ji) 
    197             IF( ht_i_b(ji) < 0.4 ) THEN 
    198                sm_i_b(ji) = 14.24 - 19.39 * ht_i_b(ji)  
    199             ELSE 
    200                sm_i_b(ji) =  7.88 - 1.59 * ht_i_b(ji) 
    201                sm_i_b(ji) = MIN( sm_i_b(ji) , zsold )   
    202             ENDIF 
    203             IF( ht_i_b(ji) > 3.06918239 ) THEN  
    204                sm_i_b(ji) = 3._wp 
    205             ENDIF 
    206             DO jk = 1, nlay_i 
    207                s_i_b(ji,jk)   = sm_i_b(ji) 
    208             END DO 
    209          END DO 
    210          ! 
    211       ENDIF ! num_sal 
     178      IF(  num_sal == 3  )   CALL lim_var_salprof1d( kideb, kiut ) 
     179 
    212180 
    213181      !------------------------------------------------------------------------------| 
    214182      ! 5) Computation of salt flux due to Bottom growth 
    215183      !------------------------------------------------------------------------------| 
    216  
    217       IF ( num_sal == 4 ) THEN 
    218          DO ji = kideb, kiut 
    219             zji = MOD( npb(ji) - 1 , jpi ) + 1 
    220             zjj =    ( npb(ji) - 1 ) / jpi + 1 
    221             fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal    )               & 
    222                &                        * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice 
    223          END DO 
    224       ELSE 
    225          DO ji = kideb, kiut 
    226             zji = MOD( npb(ji) - 1 , jpi ) + 1 
    227             zjj =    ( npb(ji) - 1 ) / jpi + 1 
    228             fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - s_i_new(ji) )               & 
    229                &                        * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice 
    230          END DO 
    231       ENDIF 
     184      ! note: s_i_new = bulk_sal in constant salinity case 
     185      DO ji = kideb, kiut 
     186         sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0._wp ) * r1_rdtice 
     187      END DO 
    232188      ! 
    233189      CALL wrk_dealloc( jpij, ze_init, zhiold, zsiold ) 

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

    r3294 r3625  
    1414   !!   lim_trp      : advection/diffusion process of sea ice 
    1515   !!---------------------------------------------------------------------- 
    16    USE phycst          ! physical constant 
    17    USE dom_oce         ! ocean domain 
    18    USE sbc_oce         ! ocean surface boundary condition 
    19    USE par_ice         ! LIM-3 parameter 
    20    USE dom_ice         ! LIM-3 domain 
    21    USE ice             ! LIM-3 variables 
    22    USE limadv          ! LIM-3 advection 
    23    USE limhdf          ! LIM-3 horizontal diffusion 
    24    USE in_out_manager  ! I/O manager 
    25    USE lbclnk          ! lateral boundary conditions -- MPP exchanges 
    26    USE lib_mpp         ! MPP library 
    27    USE wrk_nemo        ! work arrays 
    28    USE prtctl          ! Print control 
     16   USE phycst         ! physical constant 
     17   USE dom_oce        ! ocean domain 
     18   USE sbc_oce        ! ocean surface boundary condition 
     19   USE par_ice        ! ice parameter 
     20   USE dom_ice        ! ice domain 
     21   USE ice            ! ice variables 
     22   USE limadv         ! ice advection 
     23   USE limhdf         ! ice horizontal diffusion 
     24   USE in_out_manager ! I/O manager 
     25   USE lbclnk         ! lateral boundary conditions -- MPP exchanges 
     26   USE lib_mpp        ! MPP library 
     27   USE wrk_nemo       ! work arrays 
     28   USE prtctl         ! Print control 
     29   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2930 
    3031   IMPLICIT NONE 
    … …  
    4546#  include "vectopt_loop_substitute.h90" 
    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) 
    … …  
    128129         zusnit = 1.0 / REAL( initad )  
    129130         IF( zcfl > 0.5 .AND. lwp )   & 
    130             WRITE(numout,*) 'lim_trp_2 : CFL violation at day ', nday, ', cfl = ', zcfl,   & 
     131            WRITE(numout,*) 'lim_trp  : CFL violation at day ', nday, ', cfl = ', zcfl,   & 
    131132               &                        ': the ice time stepping is split in two' 
    132133 
    … …  
    174175         ELSE 
    175176            DO jk = 1, initad 
    176                CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area 
     177               CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area 
    177178                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    178                CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ow (:,:), sxopw(:,:),   & 
     179               CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ow (:,:), sxopw(:,:),   & 
    179180                  &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    180181               DO jl = 1, jpl 
    181                   CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     182                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    182183                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    183                   CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   & 
     184                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   & 
    184185                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    185                   CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     186                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    186187                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    187                   CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   & 
     188                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   & 
    188189                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    189                   CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     190                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    190191                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    191                   CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   & 
     192                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   & 
    192193                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    193194 
    194                   CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
     195                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
    195196                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    196                   CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   & 
     197                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   & 
    197198                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    198                   CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
     199                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
    199200                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    200                   CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
     201                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
    201202                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    202                   CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
     203                  CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
    203204                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    204                   CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   & 
     205                  CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   & 
    205206                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    206207                  DO layer = 1, nlay_i                                                           !--- ice heat contents --- 
    207                      CALL lim_adv_y( zusnit, v_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   &  
     208                     CALL lim_adv_y( zusnit, v_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   &  
    208209                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   & 
    209210                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 
    210                      CALL lim_adv_x( zusnit, u_ice, rone , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   &  
     211                     CALL lim_adv_x( zusnit, u_ice, rzero, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   &  
    211212                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   & 
    212213                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) ) 
    … …  
    392393 
    393394                  ! Ice salinity and age 
    394                   zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj)  , & 
    395                      zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl) 
    396                   IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &  
    397                      smv_i(ji,jj,jl) = zindic*zsal + (1.0-zindic)*0.0 
    398  
    399                   zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / &  
    400                      MAX( a_i(ji,jj,jl), epsi16 )  ), 0.0 ) * a_i(ji,jj,jl) 
    401                   oa_i (ji,jj,jl)  = zindic*zage  
     395                  zsal = MAX( MIN(  (rhoic-rhosn)/rhoic*sss_m(ji,jj) ,   & 
     396                     &              zusvoic * zs0sm(ji,jj,jl)         ) , s_i_min ) * v_i(ji,jj,jl) 
     397                  IF(  num_sal == 2  )   smv_i(ji,jj,jl) = zindic * zsal + (1.0-zindic) * 0._wp 
     398 
     399                  zage = MAX(  MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi16 ) ), 0._wp  ) * a_i(ji,jj,jl) 
     400                  oa_i (ji,jj,jl)  = zindic * zage  
    402401 
    403402                  ! Snow heat content 

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

    r3294 r3625  
    1212   !!    lim_update   : computes update of sea-ice global variables from trend terms 
    1313   !!---------------------------------------------------------------------- 
    14    USE limrhg          ! ice rheology 
    15  
    16    USE dom_oce 
    17    USE oce             ! dynamics and tracers variables 
    18    USE in_out_manager 
    19    USE sbc_oce         ! Surface boundary condition: ocean fields 
    20    USE sbc_ice         ! Surface boundary condition: ice fields 
    21    USE dom_ice 
    22    USE phycst          ! physical constants 
    23    USE ice 
    24    USE limdyn 
    25    USE limtrp 
    26    USE limthd 
    27    USE limsbc 
    28    USE limdia 
    29    USE limwri 
    30    USE limrst 
    31    USE thd_ice         ! LIM thermodynamic sea-ice variables 
    32    USE par_ice 
    33    USE limitd_th 
    34    USE limvar 
    35    USE prtctl           ! Print control 
    36    USE lbclnk           ! lateral boundary condition - MPP exchanges 
    37    USE wrk_nemo         ! work arrays 
     14   USE dom_oce        ! ocean domain 
     15   USE oce            ! dynamics and tracers variables 
     16   USE sbc_oce        ! Surface boundary condition: ocean fields 
     17   USE sbc_ice        ! Surface boundary condition: ice fields 
     18   USE phycst         ! physical constants 
     19   USE ice            ! ice variables  
     20   USE par_ice        ! ice parameters 
     21   USE thd_ice        ! ice thermodynamic variables 
     22   USE dom_ice        ! ice domain 
     23   USE limrhg         ! ice rheology 
     24   USE limdyn         ! ice dynamics 
     25   USE limtrp         ! ice transport 
     26   USE limthd         ! ice thermodynamics 
     27   USE limsbc         ! ice-oce surface boundary conditions 
     28   USE limdia         ! ice diagnostics 
     29   USE limwri         ! ice outputs 
     30   USE limrst         ! ice restart 
     31   USE limitd_th      ! ice thickness distribution (thermodynamics) 
     32   USE limvar         ! ice variables 
     33   USE prtctl         ! Print control 
     34   USE in_out_manager ! I/O manager 
     35   USE lbclnk         ! lateral boundary condition - MPP exchanges 
     36   USE wrk_nemo       ! work arrays 
     37   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3838 
    3939   IMPLICIT NONE 
    … …  
    5454#  include "vectopt_loop_substitute.h90" 
    5555   !!---------------------------------------------------------------------- 
    56    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     56   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    5757   !! $Id$ 
    5858   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    190190 
    191191                  ! is there any ice left ? 
    192                   zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )  
     192                  zindic = MAX( rzero , SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )  
    193193                  !=1 if hi > 1e-3 and 0 if not 
    194                   zdvres        = MAX(0.0,-v_i(ji,jj,jl)) !residual volume if too much ice was molten 
     194                  zdvres = MAX( 0.0 , -v_i(ji,jj,jl) ) !residual volume if too much ice was molten 
    195195                  !this quantity is positive 
    196                   v_i(ji,jj,jl) = zindic*v_i(ji,jj,jl)    !ice volume cannot be negative 
     196                  v_i(ji,jj,jl) = zindic * v_i(ji,jj,jl)    !ice volume cannot be negative 
    197197                  !correct thermodynamic ablation 
    198                   d_v_i_thd(ji,jj,jl)  = zindic  *  d_v_i_thd(ji,jj,jl) + (1.0-zindic) * (-zviold - d_v_i_trp(ji,jj,jl))  
     198                  d_v_i_thd(ji,jj,jl) = zindic  *  d_v_i_thd(ji,jj,jl) + (1.0-zindic) * (-zviold - d_v_i_trp(ji,jj,jl))  
    199199                  ! THIS IS NEW 
    200                   d_a_i_thd(ji,jj,jl)  = zindic  *  d_a_i_thd(ji,jj,jl) + &  
    201                      (1.0-zindic) * (-old_a_i(ji,jj,jl))  
     200                  d_a_i_thd(ji,jj,jl) = zindic  *  d_a_i_thd(ji,jj,jl) + (1.0-zindic) * (-old_a_i(ji,jj,jl))  
    202201 
    203202                  !residual salt flux if ice is over-molten 
    204                   fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * &  
    205                      ( rhoic * zdvres / rdt_ice ) 
    206                   !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhoic * lfus * zdvres / rdt_ice 
     203                  sfx_res(ji,jj)  = sfx_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres * r1_rdtice ) 
     204                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhoic * lfus * zdvres * r1_rdtice 
    207205 
    208206                  ! is there any snow left ? 
    209                   zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )  
    210                   zvsold        = v_s(ji,jj,jl) 
    211                   zdvres        = MAX(0.0,-v_s(ji,jj,jl)) !residual volume if too much ice was molten 
     207                  zindsn = MAX(  rzero, SIGN( rone , v_s(ji,jj,jl) - epsi10 ) )  
     208                  zvsold = v_s(ji,jj,jl) 
     209                  zdvres = MAX( 0.0 , -v_s(ji,jj,jl) )  !residual volume if too much ice was molten 
    212210                  !this quantity is positive 
    213211                  v_s(ji,jj,jl) = zindsn*v_s(ji,jj,jl)    !snow volume cannot be negative 
    … …  
    215213                  d_v_s_thd(ji,jj,jl)  = zindsn  *  d_v_s_thd(ji,jj,jl) + &  
    216214                     (1.0-zindsn) * (-zvsold - d_v_s_trp(ji,jj,jl))  
    217                   !unsure correction on salt flux.... maybe future will tell it was not that right 
    218  
    219                   !residual salt flux if snow is over-molten 
    220                   fsalt_res(ji,jj)  = fsalt_res(ji,jj) + sss_m(ji,jj) * ( rhosn * zdvres / rdt_ice ) 
    221                   !this flux will be positive if snow was over-molten 
    222                   !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhosn * lfus * zdvres / rdt_ice 
     215 
     216                  !no salt flux when snow is over-molten 
     217                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhosn * lfus * zdvres * r1_rdtice 
    223218               ENDIF 
    224219            END DO !ji 
    … …  
    229224         DO jj = 1, jpj 
    230225            DO ji = 1, jpi 
    231                IF ( ABS(fsalt_res(ji,jj)) .GT. 1.0 ) THEN  
    232                   WRITE(numout,*) ' ALERTE 1000 : residual salt flux of -> ', & 
    233                      fsalt_res(ji,jj) 
    234                   WRITE(numout,*) ' ji, jj : ', ji, jj, ' gphit, glamt : ', & 
    235                      gphit(ji,jj), glamt(ji,jj) 
     226               IF(  ABS( sfx_res(ji,jj) )  >  1._wp  ) THEN  
     227                  WRITE(numout,*) ' ALERTE 1000 : residual salt flux of -> ', sfx_res(ji,jj) 
     228                  WRITE(numout,*) ' ji, jj : ', ji, jj, ' gphit, glamt : ', gphit(ji,jj), glamt(ji,jj) 
    236229               ENDIF 
    237230            END DO 
    … …  
    277270      ENDIF 
    278271 
    279       at_i(:,:) = 0._wp 
    280       DO jl = 1, jpl 
    281          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     272      at_i(:,:) = a_i(:,:,1) 
     273      DO jl = 2, jpl 
     274         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    282275      END DO 
    283276 
    … …  
    306299      !-------------- 
    307300 
    308       IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN      ! general case 
     301      IF(  num_sal == 2  ) THEN      ! Prognostic salinity [Sice=F(z,t)] 
    309302         ! 
    310303         IF( ln_nicep ) THEN   
    … …  
    317310            WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl) 
    318311         ENDIF 
    319  
     312         ! 
    320313         smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:) + d_smv_i_trp(:,:,:) 
    321314         ! 
    … …  
    352345               .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb).GT.zhimax ) & 
    353346               .AND.( zat_i_old.LT.zacrith ) )  THEN ! new line 
    354                z_prescr_hi      = hi_max(1) / 2.0 
    355                a_i(ji,jj,1)     = v_i(ji,jj,1) / z_prescr_hi 
     347               z_prescr_hi  = hi_max(1) * 0.5_wp 
     348               a_i(ji,jj,1) = v_i(ji,jj,1) / z_prescr_hi 
    356349            ENDIF 
    357350         END DO 
    … …  
    412405      ENDIF 
    413406 
    414       at_i(:,:) = 0._wp 
    415       DO jl = 1, jpl 
    416          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     407      at_i(:,:) = a_i(:,:,1) 
     408      DO jl = 2, jpl 
     409         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    417410      END DO 
    418411 
    … …  
    452445      ENDIF 
    453446 
    454       at_i(:,:) = 0._wp 
    455       DO jl = 1, jpl 
    456          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     447      at_i(:,:) = a_i(:,:,1) 
     448      DO jl = 2, jpl 
     449         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    457450      END DO 
    458451 
    … …  
    616609            DO ji = 1, jpi 
    617610               IF ( internal_melt(ji,jj,jl) == 1 ) THEN 
    618                   v_s(ji,jj,jl)   = 0.0 
    619                   e_s(ji,jj,1,jl) = 0.0 
     611                  v_s(ji,jj,jl)   = 0._wp 
     612                  e_s(ji,jj,1,jl) = 0._wp 
    620613                  !   ! release heat 
    621                   fheat_res(ji,jj) = fheat_res(ji,jj)  & 
    622                      + ze_s * v_s(ji,jj,jl) / rdt_ice 
     614                  fheat_res(ji,jj) = fheat_res(ji,jj) + ze_s  * v_s(ji,jj,jl) * r1_rdtice 
    623615                  ! release mass 
    624                   rdmsnif(ji,jj) =  rdmsnif(ji,jj) + rhosn * v_s(ji,jj,jl) 
     616                  rdm_snw  (ji,jj) = rdm_snw  (ji,jj) + rhosn * v_s(ji,jj,jl) 
    625617               ENDIF 
    626618            END DO 
    … …  
    648640               !                ENDIF 
    649641               IF ((oa_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN 
    650                   oa_i(ji,jj,jl) = rdt_ice*numit/86400.0*a_i(ji,jj,jl) 
     642                  oa_i(ji,jj,jl) = rdt_ice * numit / 86400.0 * a_i(ji,jj,jl) 
    651643               ENDIF 
    652644               oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl) 
    … …  
    657649               !             v_s(ji,jj,jl)  = MAX( zindb * v_s(ji,jj,jl), 0.0)  
    658650               ! snow thickness cannot be smaller than 1e-6 
    659                v_s(ji,jj,jl)  = zindsn*v_s(ji,jj,jl)*zindb 
     651               v_s(ji,jj,jl)  = zindsn * v_s(ji,jj,jl) * zindb 
    660652               v_s(ji,jj,jl)  = v_s(ji,jj,jl) *  MAX( 0.0 , SIGN( 1.0 , v_s(ji,jj,jl) - epsi06 ) ) 
    661653 
    … …  
    737729      !--------------------- 
    738730 
    739       IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case 
    740  
     731      IF(  num_sal == 2  ) THEN      ! Prognostic salinity [Sice=F(z,t)] 
     732         ! 
    741733         DO jl = 1, jpl 
    742734            DO jk = 1, nlay_i 
    743735               DO jj = 1, jpj  
    744                   DO ji = 1, jpi 
    745                      ! salinity stays in bounds 
    746                      smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)), & 
    747                         0.1 * v_i(ji,jj,jl) ) 
     736                  DO ji = 1, jpi           ! salinity stays in bounds 
     737                     smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)), 0.1 * v_i(ji,jj,jl) ) 
    748738                     i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl))) 
    749                      smv_i(ji,jj,jl)  = i_ice_switch*smv_i(ji,jj,jl) + & 
    750                         0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl) 
     739                     smv_i(ji,jj,jl)  = i_ice_switch*smv_i(ji,jj,jl) + 0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl) 
    751740                  END DO ! ji 
    752741               END DO ! jj 
    753742            END DO !jk 
    754743         END DO !jl 
    755  
     744         ! 
    756745      ENDIF 
    757746 
    … …  
    796785      !----------------------------------------------------- 
    797786      zamax = amax 
    798       ! 2.13.1) individual concentrations cannot exceed zamax 
    799       !------------------------------------------------------ 
    800  
    801       at_i(:,:) = 0.0 
    802       DO jl = 1, jpl 
    803          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    804       END DO 
    805  
    806       ! 2.13.2) Total ice concentration cannot exceed zamax 
    807       !---------------------------------------------------- 
     787      ! 2.13.1) total (and thus individual) concentrations cannot exceed zamax 
     788      !----------------------------------------------------------------------- 
     789 
    808790      at_i(:,:) = a_i(:,:,1) 
    809791      DO jl = 2, jpl 
    810          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     792         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    811793      END DO 
    812794 
    … …  
    815797 
    816798            ! 0) Excessive area ? 
    817             z_da_ex =  MAX( at_i(ji,jj) - zamax , 0.0 )         
     799            z_da_ex =  MAX( at_i(ji,jj) - zamax , 0._wp )         
    818800 
    819801            ! 1) Count the number of existing categories 
    820802            DO jl  = 1, jpl 
     803!!cr : comment the second line of zindb definition, and use epsi04 in the 1st one 
    821804               zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi03 ) )  
    822805               zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) ) )  
    … …  
    839822      at_i(:,:) = a_i(:,:,1) 
    840823      DO jl = 2, jpl 
    841          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     824         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    842825      END DO 
    843826 
    … …  
    894877      at_i(:,:) = a_i(:,:,1) 
    895878      DO jl = 2, jpl 
    896          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     879         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    897880      END DO 
    898881 
    … …  
    902885      ! Ice drift 
    903886      !------------ 
     887!!gm  BUG ? I don't understand this : it may have a wrong impact on the ice edge advection 
     888!!gm  and any way there is much faster way to code that... 
    904889      DO jj = 2, jpjm1 
    905890         DO ji = fs_2, fs_jpim1 
    … …  
    913898      END DO 
    914899      !mask velocities 
     900!!gm BUG ?  here the mask are the one of the beginning of the time step, no?  
     901!!gm        whereas at this level they should have been updated... To be checked  
    915902      u_ice(:,:) = u_ice(:,:) * tmu(:,:) 
    916903      v_ice(:,:) = v_ice(:,:) * tmv(:,:) 
    … …  
    10211008         CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ') 
    10221009         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ') 
    1023          CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update : fhbri : ', tab2d_2=fheat_rpo , clinfo2= ' fheat_rpo : ') 
     1010         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update : fhbri : ', tab2d_2=fheat_mec , clinfo2= ' fheat_mec : ') 
    10241011 
    10251012         CALL prt_ctl_info(' ') 

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

    r3294 r3625  
    4343   !!   lim_var_bv        : 
    4444   !!---------------------------------------------------------------------- 
    45    USE par_oce          ! ocean parameters 
    46    USE phycst           ! physical constants (ocean directory)  
    47    USE sbc_oce          ! Surface boundary condition: ocean fields 
    48    USE ice              ! LIM variables 
    49    USE par_ice          ! LIM parameters 
    50    USE dom_ice          ! LIM domain 
    51    USE thd_ice          ! LIM thermodynamics 
    52    USE in_out_manager   ! I/O manager 
    53    USE lib_mpp          ! MPP library 
    54    USE wrk_nemo         ! work arrays 
     45   USE par_oce        ! ocean parameters 
     46   USE phycst         ! physical constants (ocean directory)  
     47   USE sbc_oce        ! Surface boundary condition: ocean fields 
     48   USE ice            ! ice variables 
     49   USE par_ice        ! ice parameters 
     50   USE thd_ice        ! ice variables (thermodynamics) 
     51   USE dom_ice        ! ice domain 
     52   USE in_out_manager ! I/O manager 
     53   USE lib_mpp        ! MPP library 
     54   USE wrk_nemo       ! work arrays 
     55   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    5556 
    5657   IMPLICIT NONE 
    … …  
    7374 
    7475   !!---------------------------------------------------------------------- 
    75    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     76   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    7677   !! $Id$ 
    7778   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    182183      END DO 
    183184 
    184       IF(  num_sal == 2  .OR.  num_sal == 4  )THEN 
     185      IF(  num_sal == 2  )THEN 
    185186         DO jl = 1, jpl 
    186187            DO jj = 1, jpj 
    … …  
    309310      ! Vertically constant, constant in time 
    310311      !--------------------------------------- 
    311       IF( num_sal == 1 )   s_i(:,:,:,:) = bulk_sal 
     312      IF(  num_sal == 1 )   s_i(:,:,:,:) = bulk_sal 
    312313 
    313314      !----------------------------------- 
    314315      ! Salinity profile, varying in time 
    315316      !----------------------------------- 
    316  
    317       IF(   num_sal == 2  .OR.   num_sal == 4   ) THEN 
     317      IF(  num_sal == 2  ) THEN 
    318318         ! 
    319319         DO jk = 1, nlay_i 
    … …  
    331331         dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 )       ! Weighting factor between zs_zero and zs_inf 
    332332         dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 ) 
    333  
     333         ! 
    334334         zalpha(:,:,:) = 0._wp 
    335335         DO jl = 1, jpl 
    … …  
    347347            END DO 
    348348         END DO 
    349  
     349         ! 
    350350         dummy_fac = 1._wp / nlay_i                   ! Computation of the profile 
    351351         DO jl = 1, jpl 
    … …  
    361361            END DO ! jk 
    362362         END DO ! jl 
    363  
     363         ! 
    364364      ENDIF ! num_sal 
    365365 
    … …  
    368368      !------------------------------------------------------- 
    369369 
    370       IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 
     370      IF(  num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 
    371371         ! 
    372372         sm_i(:,:,:) = 2.30_wp 
    … …  
    380380            END DO 
    381381         END DO 
    382  
     382         ! 
    383383      ENDIF ! num_sal 
    384384      ! 
    … …  
    447447      !------------------------------------------------------ 
    448448 
    449       IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN 
     449      IF(  num_sal == 2  ) THEN 
    450450         ! 
    451451         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero 

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

    r3294 r3625  
    2525   USE wrk_nemo        ! work arrays 
    2626   USE par_ice 
     27   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2728 
    2829   IMPLICIT NONE 
    … …  
    5152   REAL(wp)  ::   zone   = 1._wp       
    5253   !!---------------------------------------------------------------------- 
    53    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     54   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    5455   !! $Id$ 
    5556   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    245246            zcmo(ji,jj,25) = et_i(ji,jj) 
    246247            zcmo(ji,jj,26) = et_s(ji,jj) 
    247             zcmo(ji,jj,28) = fsbri(ji,jj) 
    248             zcmo(ji,jj,29) = fseqv(ji,jj) 
     248            zcmo(ji,jj,28) = sfx_bri(ji,jj) 
     249            zcmo(ji,jj,29) = sfx_thd(ji,jj) 
    249250 
    250251            zcmo(ji,jj,30) = bv_i(ji,jj) 

  • branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90

    r2715 r3625  
    111111         zcmo(ji,jj,13) = qns(ji,jj) 
    112112         ! See thersf for the coefficient 
    113          zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     113         zcmo(ji,jj,14) = - sfx (ji,jj) * rday      ! converted in Kg/m2/day = mm/day 
    114114         zcmo(ji,jj,15) = utau_ice(ji,jj) 
    115115         zcmo(ji,jj,16) = vtau_ice(ji,jj) 
    … …  
    154154               rcmoy(ji,jj,13) = qns(ji,jj) 
    155155               ! See thersf for the coefficient 
    156                rcmoy(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     156               rcmoy(ji,jj,14) = - sfx (ji,jj) * rday      ! converted in mm/day 
    157157               rcmoy(ji,jj,15) = utau_ice(ji,jj) 
    158158               rcmoy(ji,jj,16) = vtau_ice(ji,jj) 

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

    r2715 r3625  
    88   USE par_ice        ! LIM-3 parameters 
    99   USE in_out_manager ! I/O manager 
    10    USE lib_mpp         ! MPP library 
     10   USE lib_mpp        ! MPP library 
    1111 
    1212   IMPLICIT NONE 
    … …  
    6666   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_b        !: <==> the 2D  frld 
    6767   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fbif_1d       !: <==> the 2D  fbif 
    68    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdmicif_1d    !: <==> the 2D  rdmicif 
    69    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdmsnif_1d    !: <==> the 2D  rdmsnif 
     68   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdm_ice_1d    !: <==> the 2D  rdm_ice 
     69   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdm_snw_1d    !: <==> the 2D  rdm_snw 
    7070   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlbbq_1d      !: <==> the 2D  qlbsbq 
    7171   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dmgwi_1d      !: <==> the 2D  dmgwi 
    … …  
    8383   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice 
    8484   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   old_ht_i_b    !: Ice thickness at the beginnning of the time step [m] 
    85     REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::  old_ht_s_b    !: Snow thickness at the beginning of the time step [m] 
    86    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fsbri_1d      !: Salt flux due to brine drainage 
     85   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   old_ht_s_b    !: Snow thickness at the beginning of the time step [m] 
     86   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_bri_1d    !: <==> the 2D sfx_bri 
    8787   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhbri_1d      !: Heat flux due to brine drainage 
    88    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fseqv_1d      !: Equivalent Salt flux due to ice growth/decay 
     88   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_thd_1d    !: <==> the 2D sfx_thd 
    8989   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dsm_i_fl_1d   !: Ice salinity variations due to flushing 
    9090   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dsm_i_gd_1d   !: Ice salinity variations due to gravity drainage 
    … …  
    138138 
    139139   !!---------------------------------------------------------------------- 
    140    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     140   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
    141141   !! $Id$ 
    142142   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    160160      ! 
    161161      ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_b     (jpij) ,     & 
    162          &      fbif_1d    (jpij) , rdmicif_1d (jpij) , rdmsnif_1d (jpij) ,     & 
     162         &      fbif_1d    (jpij) , rdm_ice_1d (jpij) , rdm_snw_1d (jpij) ,     & 
    163163         &      qlbbq_1d   (jpij) , dmgwi_1d   (jpij) , dvsbq_1d   (jpij) ,     & 
    164164         &      dvbbq_1d   (jpij) , dvlbq_1d   (jpij) , dvnbq_1d   (jpij) ,     & 
    … …  
    166166         &      tatm_ice_1d(jpij) , fsup       (jpij) , focea      (jpij) ,     &    
    167167         &      i0         (jpij) , old_ht_i_b (jpij) , old_ht_s_b (jpij) ,     &   
    168          &      fsbri_1d   (jpij) , fhbri_1d   (jpij) , fseqv_1d  (jpij) ,     & 
     168         &      sfx_bri_1d (jpij) , fhbri_1d   (jpij) , sfx_thd_1d (jpij) ,     & 
    169169         &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,     &      
    170170         &      dsm_i_si_1d(jpij) , hicol_b    (jpij)                     , STAT=ierr(2) ) 
Note: See TracChangeset for help on using the changeset viewer.