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

Changeset 4099


Ignore:
Timestamp:
2013-10-22T14:07:21+02:00 (10 years ago)
Author:
clem
Message:
 
Location:
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM
Files:
2 added
2 deleted
17 edited

Legend:

Unmodified
Added
Removed

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/cfg.txt

    r3905 r4099  
    1 AMM12 OPA_SRC 
    21GYRE OPA_SRC 
    32GYRE_BFM OPA_SRC TOP_SRC 
    43GYRE_PISCES OPA_SRC TOP_SRC 
    5 ORCA2_LIM3 OPA_SRC LIM_SRC_3 
    64ORCA2_SAS_LIM OPA_SRC SAS_SRC LIM_SRC_2 NST_SRC 
    75ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    8 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    96ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
    107ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
     8AMM12 OPA_SRC 
     9ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
     10ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r4045 r4099  
    395395   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out = "restart_ice"      !: suffix of ice restart name (output) 
    396396   LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F) 
    397    LOGICAL               , PUBLIC ::   ln_nicep      = .TRUE.             !: flag for sea-ice points output (T) or not (F) 
     397   LOGICAL               , PUBLIC ::   ln_nicep      = .FALSE.             !: flag for sea-ice points output (T) or not (F) 
    398398   REAL(wp)              , PUBLIC ::   cai           = 1.40e-3            !: atmospheric drag over sea ice 
    399399   REAL(wp)              , PUBLIC ::   cao           = 1.00e-3            !: atmospheric drag over ocean 
    … …  
    404404   !!-------------------------------------------------------------------------- 
    405405   !! Check if everything down here is necessary 
    406    LOGICAL , PUBLIC                                      ::   ln_limdiahsb  = .TRUE. !: flag for ice diag (T) or not (F) 
     406   LOGICAL , PUBLIC                                      ::   ln_limdiahsb  = .FALSE. !: flag for ice diag (T) or not (F) 
     407   LOGICAL , PUBLIC                                      ::   ln_limdiaout  = .FALSE. !: flag for ice diag (T) or not (F) 
    407408   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   v_newice   !: volume of ice formed in the leads 
    408409   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dv_dt_thd  !: thermodynamic growth rates  

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90

    r4045 r4099  
    126126      !! ** input   :   Namelist namicerun 
    127127      !!------------------------------------------------------------------- 
    128       NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb 
     128      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb, ln_limdiaout 
    129129      !!------------------------------------------------------------------- 
    130130      !                     
    … …  
    147147         WRITE(numout,*) '   Several ice points in the ice or not in ocean.output    = ', ln_nicep 
    148148         WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb 
     149         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout 
    149150      ENDIF 
    150151      ! 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90

    r4072 r4099  
    154154                                                           !============================ 
    155155         ! Check if tests have passed (i.e. volume conservation...) 
    156          IF ( ztests /= 4 ) THEN 
    157             WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
    158             WRITE(numout,*) ' !! ALERT categories distribution !!' 
    159             WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
    160             WRITE(numout,*) ' *** ztests is not equal to 4 ' 
    161             WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
    162             WRITE(numout,*) 'i_fill=',i_fill 
    163             WRITE(numout,*) 'zai(ji)=',zai(ji) 
    164             WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:) 
    165          ENDIF 
     156         !IF ( ztests /= 4 ) THEN 
     157         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
     158         !   WRITE(numout,*) ' !! ALERT categories distribution !!' 
     159         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
     160         !   WRITE(numout,*) ' *** ztests is not equal to 4 ' 
     161         !   WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
     162         !   WRITE(numout,*) 'i_fill=',i_fill 
     163         !   WRITE(numout,*) 'zai(ji)=',zai(ji) 
     164         !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:) 
     165         !ENDIF 
    166166             
    167167      END DO ! i loop 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_2D.F90

    r4072 r4099  
    155155                                                              !============================ 
    156156            ! Check if tests have passed (i.e. volume conservation...) 
    157             IF ( ztests .NE. 4 ) THEN 
    158                WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
    159                WRITE(numout,*) ' !! ALERT categories distribution !!' 
    160                WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
    161                WRITE(numout,*) ' *** ztests is not equal to 4 ' 
    162                WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
    163             ENDIF 
     157            !IF ( ztests .NE. 4 ) THEN 
     158            !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
     159            !   WRITE(numout,*) ' !! ALERT categories distribution !!' 
     160            !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
     161            !   WRITE(numout,*) ' *** ztests is not equal to 4 ' 
     162            !   WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
     163            !ENDIF 
    164164 
    165165         END DO ! i loop 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r4045 r4099  
    1919   USE dom_oce          ! ocean domain 
    2020   USE sbc_oce          ! Surface boundary condition: ocean fields 
     21   USE sbc_ice          ! Surface boundary condition: ice fields 
    2122   USE eosbn2           ! equation of state 
    2223   USE ice              ! sea-ice variables 
    … …  
    2728   USE lbclnk           ! lateral boundary condition - MPP exchanges 
    2829   USE lib_mpp          ! MPP library 
     30   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2931   USE wrk_nemo         ! work arrays 
    30    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3132 
    3233   IMPLICIT NONE 
    … …  
    3536   PUBLIC   lim_istate      ! routine called by lim_init.F90 
    3637 
    37    !                                  !!** init namelist (namiceini) ** 
    38    REAL(wp) ::   ttest    = 2.0_wp     ! threshold water temperature for initial sea ice 
    39    REAL(wp) ::   hninn    = 0.5_wp     ! initial snow thickness in the north 
    40    REAL(wp) ::   hginn_u  = 2.5_wp     ! initial ice thickness in the north 
    41    REAL(wp) ::   aginn_u  = 0.7_wp     ! initial leads area in the north 
    42    REAL(wp) ::   hginn_d  = 5.0_wp     ! initial ice thickness in the north 
    43    REAL(wp) ::   aginn_d  = 0.25_wp    ! initial leads area in the north 
    44    REAL(wp) ::   hnins    = 0.1_wp     ! initial snow thickness in the south 
    45    REAL(wp) ::   hgins_u  = 1.0_wp     ! initial ice thickness in the south 
    46    REAL(wp) ::   agins_u  = 0.7_wp     ! initial leads area in the south 
    47    REAL(wp) ::   hgins_d  = 2.0_wp     ! initial ice thickness in the south 
    48    REAL(wp) ::   agins_d  = 0.2_wp     ! initial leads area in the south 
    49    REAL(wp) ::   sinn     = 6.301_wp   ! initial salinity  
    50    REAL(wp) ::   sins     = 6.301_wp   ! 
    51  
    52    !!---------------------------------------------------------------------- 
    53    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     38   !! * Module variables 
     39   REAL(wp) ::             & !!! ** init namelist (namiceini) ** 
     40      ttest    = 2.0  ,    &  ! threshold water temperature for initial sea ice 
     41      hninn    = 0.5  ,    &  ! initial snow thickness in the north 
     42      hnins    = 0.1  ,    &  ! initial snow thickness in the south 
     43      hginn    = 2.5  ,    &  ! initial ice thickness in the north 
     44      hgins    = 1.0  ,    &  ! initial ice thickness in the south 
     45      aginn    = 0.7  ,    &  ! initial leads area in the north 
     46      agins    = 0.7  ,    &  ! initial leads area in the south 
     47      sinn     = 6.301 ,   &  ! initial salinity  
     48      sins     = 6.301 
     49 
     50   !!---------------------------------------------------------------------- 
     51   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008) 
    5452   !! $Id$ 
    55    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    56    !!---------------------------------------------------------------------- 
     53   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     54   !!---------------------------------------------------------------------- 
     55 
    5756CONTAINS 
    5857 
    … …  
    6362      !! ** Purpose :   defined the sea-ice initial state 
    6463      !! 
    65       !! ** Method  :   restart from a state defined in a binary file 
    66       !!                or from arbitrary sea-ice conditions 
    67       !!------------------------------------------------------------------- 
    68       INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices 
    69       REAL(wp) ::   zeps6, zeps, ztmelts, epsi06   ! local scalars 
    70       REAL(wp) ::   zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs  
    71       REAL(wp), POINTER, DIMENSION(:)   ::   zgfactorn, zhin  
    72       REAL(wp), POINTER, DIMENSION(:)   ::   zgfactors, zhis 
    73       REAL(wp), POINTER, DIMENSION(:,:) ::   zidto      ! ice indicator 
    74       !-------------------------------------------------------------------- 
    75  
    76       CALL wrk_alloc( jpm, zgfactorn, zgfactors, zhin, zhis ) 
     64      !! ** Method  :    
     65      !!                This routine will put some ice where ocean 
     66      !!                is at the freezing point, then fill in ice  
     67      !!                state variables using prescribed initial  
     68      !!                values in the namelist             
     69      !! 
     70      !! ** Steps   :    
     71      !!                1) Read namelist 
     72      !!                2) Basal temperature; ice and hemisphere masks 
     73      !!                3) Fill in the ice thickness distribution using gaussian 
     74      !!                4) Fill in space-dependent arrays for state variables 
     75      !!                5) Diagnostic arrays 
     76      !!                6) Lateral boundary conditions 
     77      !! 
     78      !! History : 
     79      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code 
     80      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats 
     81      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats 
     82      !!-------------------------------------------------------------------- 
     83 
     84      !! * Local variables 
     85      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices 
     86      REAL(wp)   :: epsi06, epsi20, ztmelts 
     87      INTEGER    :: i_hemis, i_fill, jl0   
     88      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 
     89      REAL(wp), POINTER, DIMENSION(:)     :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini 
     90      REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, za_i_ini, zv_i_ini 
     91      REAL(wp), POINTER, DIMENSION(:,:)   :: zidto    ! ice indicator 
     92      INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index 
     93      !-------------------------------------------------------------------- 
     94 
    7795      CALL wrk_alloc( jpi, jpj, zidto ) 
    78  
    79       !-------------------------------------------------------------------- 
    80       ! 1) Preliminary things  
    81       !-------------------------------------------------------------------- 
    82       epsi06 = 1.e-6_wp 
     96      CALL wrk_alloc( jpi, jpj, zhemis ) 
     97      CALL wrk_alloc( jpl,   2, zht_i_ini,  za_i_ini,  zv_i_ini ) 
     98      CALL wrk_alloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 
     99 
     100      epsi06   = 1.0e-6 
     101      epsi20   = 1.0e-20 
     102      IF(lwp) WRITE(numout,*) 
     103      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' 
     104      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 
     105 
     106      !-------------------------------------------------------------------- 
     107      ! 1) Read namelist 
     108      !-------------------------------------------------------------------- 
    83109 
    84110      CALL lim_istate_init     !  reading the initials parameters of the ice 
    … …  
    89115 
    90116      !-------------------------------------------------------------------- 
    91       ! 2) Ice initialization (hi,hs,frld,t_su,sm_i,t_i,t_s)              |  
    92       !-------------------------------------------------------------------- 
    93  
    94       IF(lwp) WRITE(numout,*) 
    95       IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' 
    96       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 
    97  
     117      ! 2) Basal temperature, ice mask and hemispheric index 
     118      !-------------------------------------------------------------------- 
     119 
     120      ! Basal temperature is set to the freezing point of seawater in Celsius 
    98121      t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
    99122 
    100123      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest 
    101124         DO ji = 1, jpi 
    102             IF( tsn(ji,jj,1,jp_tem)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0.e0      ! no ice 
    103             ELSE                                                     ;   zidto(ji,jj) = 1.e0      !    ice 
     125            IF( tsn(ji,jj,1,jp_tem)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0._wp      ! no ice 
     126            ELSE                                                     ;   zidto(ji,jj) = 1._wp      !    ice 
    104127            ENDIF 
    105128         END DO 
    106129      END DO 
    107130 
    108       t_bo(:,:) = t_bo(:,:) + rt0                          ! t_bo converted from Celsius to Kelvin (rt0 over land) 
    109  
    110       ! constants for heat contents 
    111       zeps   = 1.e-20_wp 
    112       zeps6  = 1.e-06_wp 
    113  
    114       ! zgfactor for initial ice distribution 
    115       zgfactorn(:) = 0._wp 
    116       zgfactors(:) = 0._wp 
    117  
    118       ! first ice type 
    119       DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
    120          zhin (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 
    121          zgfactorn(1) = zgfactorn(1) + exp(-(zhin(1)-hginn_u)*(zhin(1)-hginn_u) * 0.5_wp ) 
    122          zhis (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 
    123          zgfactors(1) = zgfactors(1) + exp(-(zhis(1)-hgins_u)*(zhis(1)-hgins_u) * 0.5_wp ) 
    124       END DO ! jl 
    125       zgfactorn(1) = aginn_u / zgfactorn(1) 
    126       zgfactors(1) = agins_u / zgfactors(1) 
    127  
    128       ! ------------- 
    129       ! new distribution, polynom of second order, conserving area and volume 
    130       zh1 = 0._wp 
    131       zh2 = 0._wp 
    132       zh3 = 0._wp 
    133       DO jl = 1, jpl 
    134          zh = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 
    135          zh1 = zh1 + zh 
    136          zh2 = zh2 + zh * zh 
    137          zh3 = zh3 + zh * zh * zh 
    138       END DO 
    139       IF(lwp) WRITE(numout,*) ' zh1 : ', zh1 
    140       IF(lwp) WRITE(numout,*) ' zh2 : ', zh2 
    141       IF(lwp) WRITE(numout,*) ' zh3 : ', zh3 
    142  
    143       zvol = aginn_u * hginn_u 
    144       zare = aginn_u 
    145       IF( jpl >= 2 ) THEN 
    146          zbn = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3) 
    147          zan = ( zare - zbn*zh1 ) / zh2 
    148       ENDIF 
    149  
    150       IF(lwp) WRITE(numout,*) ' zvol: ', zvol 
    151       IF(lwp) WRITE(numout,*) ' zare: ', zare 
    152       IF(lwp) WRITE(numout,*) ' zbn : ', zbn  
    153       IF(lwp) WRITE(numout,*) ' zan : ', zan  
    154  
    155       zvol = agins_u * hgins_u 
    156       zare = agins_u 
    157       IF( jpl >= 2 ) THEN 
    158          zbs = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3) 
    159          zas = ( zare - zbs*zh1 ) / zh2 
    160       ENDIF 
    161  
    162       IF(lwp) WRITE(numout,*) ' zvol: ', zvol 
    163       IF(lwp) WRITE(numout,*) ' zare: ', zare 
    164       IF(lwp) WRITE(numout,*) ' zbn : ', zbn  
    165       IF(lwp) WRITE(numout,*) ' zan : ', zan  
    166  
    167       !end of new lines 
    168       ! ------------- 
    169 !!! 
    170       ! retour a LIMA_MEC 
    171       !     ! second ice type 
    172       !     zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
    173       !     hi_max(ice_cat_bounds(2,1)-1) = 0.0 
    174  
    175       !     ! here to change !!!! 
    176       !     jm = 2 
    177       !     DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    178       !        zhin (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    179       !        zhin (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + & 
    180       !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0 
    181       !        zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0) 
    182       !        zhis (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    183       !        zhis (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + & 
    184       !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0 
    185       !        zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0) 
    186       !     END DO ! jl 
    187       !     zgfactorn(2) = aginn_d / zgfactorn(2) 
    188       !     zgfactors(2) = agins_d / zgfactors(2) 
    189       !     hi_max(ice_cat_bounds(2,1)-1) = zdummy 
    190       ! END retour a LIMA_MEC 
    191 !!! 
    192  
    193 !!gm  optimisation :  loop over the ice categories inside the ji, jj loop !!! 
    194  
     131      t_bo(:,:) = t_bo(:,:) + rt0                          ! conversion to Kelvin 
     132 
     133      ! Hemispheric index 
     134      ! MV 2011 new initialization 
    195135      DO jj = 1, jpj 
    196136         DO ji = 1, jpi 
    197  
    198             !--- Northern hemisphere 
    199             !---------------------------------------------------------------- 
    200137            IF( fcor(ji,jj) >= 0._wp ) THEN     
    201  
    202                !----------------------- 
    203                ! Ice area / thickness 
    204                !----------------------- 
    205  
    206                IF ( jpl .EQ. 1) THEN ! one category 
    207  
    208                   DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 
    209                      a_i(ji,jj,jl)    = zidto(ji,jj) * aginn_u 
    210                      ht_i(ji,jj,jl)   = zidto(ji,jj) * hginn_u 
    211                      v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    212                   END DO 
    213  
    214                ELSE ! several categories 
    215  
    216                   DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 
    217                      zhin(1)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    218                      a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* &  
    219                         (zhin(1)-hginn_u)/2.0) , epsi06) 
    220                      ! new line 
    221                      a_i(ji,jj,jl)    = zidto(ji,jj) * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) ) 
    222                      ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(1)  
    223                      v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    224                   END DO 
    225  
    226                ENDIF 
    227  
    228  
    229 !!! 
    230                ! retour a LIMA_MEC 
    231                !              !ridged ice 
    232                !              zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
    233                !              hi_max(ice_cat_bounds(2,1)-1) = 0.0 
    234                !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories 
    235                !                 zhin(2)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    236                !                 a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* & 
    237                !                                    (zhin(2)-hginn_d)/2.0) , epsi06) 
    238                !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(2)  
    239                !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    240                !              END DO 
    241                !              hi_max(ice_cat_bounds(2,1)-1) = zdummy 
    242  
    243                !              !rafted ice 
    244                !              jl = 6 
    245                !              a_i(ji,jj,jl)       = 0.0 
    246                !              ht_i(ji,jj,jl)      = 0.0 
    247                !              v_i(ji,jj,jl)       = 0.0 
    248                ! END retour a LIMA_MEC 
    249 !!! 
    250  
    251                DO jl = 1, jpl 
    252  
    253                   !------------- 
    254                   ! Snow depth 
    255                   !------------- 
    256                   ht_s(ji,jj,jl)   = zidto(ji,jj) * hninn 
    257                   v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 
    258  
    259                   !--------------- 
    260                   ! Ice salinity 
    261                   !--------------- 
    262                   sm_i(ji,jj,jl)   = zidto(ji,jj) * sinn  + ( 1.0 - zidto(ji,jj) ) * 0.1 
    263                   smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 
    264  
    265                   !---------- 
    266                   ! Ice age 
    267                   !---------- 
    268                   o_i(ji,jj,jl)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) ) 
    269                   oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl) 
    270  
    271                   !------------------------------ 
    272                   ! Sea ice surface temperature 
    273                   !------------------------------ 
    274  
    275                   t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 
    276  
    277                   !------------------------------------ 
    278                   ! Snow temperature and heat content 
    279                   !------------------------------------ 
    280  
    281                   DO jk = 1, nlay_s 
    282                      t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
    283                      ! Snow energy of melting 
    284                      e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
    285                      ! Change dimensions 
    286                      e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
    287                      ! Multiply by volume, so that heat content in 10^9 Joules 
    288                      e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 
    289                         v_s(ji,jj,jl)  / nlay_s 
    290                   END DO !jk 
    291  
    292                   !----------------------------------------------- 
    293                   ! Ice salinities, temperature and heat content  
    294                   !----------------------------------------------- 
    295  
    296                   DO jk = 1, nlay_i 
    297                      t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt  
    298                      s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1 
    299                      ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
    300  
    301                      ! heat content per unit volume 
    302                      e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 
    303                         (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    304                         +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 
    305                         - rcp      * ( ztmelts - rtt ) & 
    306                         ) 
    307  
    308                      ! Correct dimensions to avoid big values 
    309                      e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    310  
    311                      ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    312                      e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &  
    313                         area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 
    314                         nlay_i 
    315                   END DO ! jk 
    316  
    317                END DO ! jl  
    318  
    319             ELSE ! on fcor  
    320  
    321                !--- Southern hemisphere 
    322                !---------------------------------------------------------------- 
    323  
    324                !----------------------- 
    325                ! Ice area / thickness 
    326                !----------------------- 
    327  
    328                IF ( jpl .EQ. 1) THEN ! one category 
    329  
    330                   DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 
    331                      a_i(ji,jj,jl)    = zidto(ji,jj) * agins_u 
    332                      ht_i(ji,jj,jl)   = zidto(ji,jj) * hgins_u 
    333                      v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    334                   END DO 
    335  
    336                ELSE ! several categories 
    337  
    338                   !level ice 
    339                   DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories 
    340  
    341                      zhis(1)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    342                      a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * &  
    343                         (zhis(1)-hgins_u)/2.0) , epsi06 ) 
    344                      ! new line square distribution volume conserving 
    345                      a_i(ji,jj,jl)    = zidto(ji,jj) * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) ) 
    346                      ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(1)  
    347                      v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    348  
    349                   END DO ! jl 
    350  
    351                ENDIF 
    352  
    353 !!! 
    354                ! retour a LIMA_MEC 
    355                !              !ridged ice 
    356                !              zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
    357                !              hi_max(ice_cat_bounds(2,1)-1) = 0.0 
    358                !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories 
    359                !                 zhis(2)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    360                !                 a_i(ji,jj,jl) = zidto(ji,jj)*MAX( zgfactors(2)   & 
    361                !                    &          * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 ) 
    362                !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(2)  
    363                !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    364                !              END DO 
    365                !              hi_max(ice_cat_bounds(2,1)-1) = zdummy 
    366  
    367                !              !rafted ice 
    368                !              jl = 6 
    369                !              a_i(ji,jj,jl)       = 0.0 
    370                !              ht_i(ji,jj,jl)      = 0.0 
    371                !              v_i(ji,jj,jl)       = 0.0 
    372                ! END retour a LIMA_MEC 
    373 !!! 
    374  
    375                DO jl = 1, jpl !over thickness categories 
    376  
    377                   !--------------- 
    378                   ! Snow depth 
    379                   !--------------- 
    380  
    381                   ht_s(ji,jj,jl)   = zidto(ji,jj) * hnins 
    382                   v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 
    383  
    384                   !--------------- 
    385                   ! Ice salinity 
    386                   !--------------- 
    387  
    388                   sm_i(ji,jj,jl)   = zidto(ji,jj) * sins  + ( 1.0 - zidto(ji,jj) ) * 0.1 
    389                   smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 
    390  
    391                   !---------- 
    392                   ! Ice age 
    393                   !---------- 
    394  
    395                   o_i(ji,jj,jl)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) ) 
    396                   oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl) 
    397  
    398                   !------------------------------ 
    399                   ! Sea ice surface temperature 
    400                   !------------------------------ 
    401  
    402                   t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 
    403  
    404                   !---------------------------------- 
    405                   ! Snow temperature / heat content 
    406                   !---------------------------------- 
    407  
    408                   DO jk = 1, nlay_s 
    409                      t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
    410                      ! Snow energy of melting 
    411                      e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
    412                      ! Change dimensions 
    413                      e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
    414                      ! Multiply by volume, so that heat content in 10^9 Joules 
    415                      e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 
    416                         v_s(ji,jj,jl)  / nlay_s 
    417                   END DO 
    418  
    419                   !--------------------------------------------- 
    420                   ! Ice temperature, salinity and heat content 
    421                   !--------------------------------------------- 
    422  
    423                   DO jk = 1, nlay_i 
    424                      t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt  
    425                      s_i(ji,jj,jk,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1 
    426                      ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
    427  
    428                      ! heat content per unit volume 
    429                      e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 
    430                         (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    431                         +   lfus  * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 
    432                         - rcp      * ( ztmelts - rtt ) & 
    433                         ) 
    434  
    435                      ! Correct dimensions to avoid big values 
    436                      e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    437  
    438                      ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    439                      e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &  
    440                         area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 
    441                         nlay_i 
    442                   END DO !jk 
    443  
    444                END DO ! jl  
    445  
    446             ENDIF ! on fcor 
    447  
     138               zhemis(ji,jj) = 1 ! Northern hemisphere 
     139            ELSE 
     140               zhemis(ji,jj) = 2 ! Southern hemisphere 
     141            ENDIF 
    448142         END DO 
    449143      END DO 
    450  
    451       !-------------------------------------------------------------------- 
    452       ! 3) Global ice variables for output diagnostics                    |  
    453       !-------------------------------------------------------------------- 
    454  
    455       fsbbq (:,:)     = 0.e0 
    456       u_ice (:,:)     = 0.e0 
    457       v_ice (:,:)     = 0.e0 
    458       stress1_i(:,:)  = 0.0 
    459       stress2_i(:,:)  = 0.0 
    460       stress12_i(:,:) = 0.0 
    461  
    462       !-------------------------------------------------------------------- 
    463       ! 4) Moments for advection 
    464       !-------------------------------------------------------------------- 
    465  
    466       sxopw (:,:) = 0.e0  
    467       syopw (:,:) = 0.e0  
    468       sxxopw(:,:) = 0.e0  
    469       syyopw(:,:) = 0.e0  
    470       sxyopw(:,:) = 0.e0 
    471  
    472       sxice (:,:,:)  = 0.e0   ;   sxsn (:,:,:)  = 0.e0   ;   sxa  (:,:,:)  = 0.e0 
    473       syice (:,:,:)  = 0.e0   ;   sysn (:,:,:)  = 0.e0   ;   sya  (:,:,:)  = 0.e0 
    474       sxxice(:,:,:)  = 0.e0   ;   sxxsn(:,:,:)  = 0.e0   ;   sxxa (:,:,:)  = 0.e0 
    475       syyice(:,:,:)  = 0.e0   ;   syysn(:,:,:)  = 0.e0   ;   syya (:,:,:)  = 0.e0 
    476       sxyice(:,:,:)  = 0.e0   ;   sxysn(:,:,:)  = 0.e0   ;   sxya (:,:,:)  = 0.e0 
    477  
    478       sxc0  (:,:,:)  = 0.e0   ;   sxe  (:,:,:,:)= 0.e0    
    479       syc0  (:,:,:)  = 0.e0   ;   sye  (:,:,:,:)= 0.e0    
    480       sxxc0 (:,:,:)  = 0.e0   ;   sxxe (:,:,:,:)= 0.e0    
    481       syyc0 (:,:,:)  = 0.e0   ;   syye (:,:,:,:)= 0.e0    
    482       sxyc0 (:,:,:)  = 0.e0   ;   sxye (:,:,:,:)= 0.e0    
    483  
    484       sxsal  (:,:,:)  = 0.e0 
    485       sysal  (:,:,:)  = 0.e0 
    486       sxxsal (:,:,:)  = 0.e0 
    487       syysal (:,:,:)  = 0.e0 
    488       sxysal (:,:,:)  = 0.e0 
    489  
    490       !-------------------------------------------------------------------- 
    491       ! 5) Lateral boundary conditions                                    |  
     144      ! END MV 2011 new initialization 
     145 
     146      !-------------------------------------------------------------------- 
     147      ! 3) Initialization of sea ice state variables 
     148      !-------------------------------------------------------------------- 
     149 
     150      !----------------------------- 
     151      ! 3.1) Hemisphere-dependent arrays 
     152      !----------------------------- 
     153      ! assign initial thickness, concentration, snow depth and salinity to 
     154      ! an hemisphere-dependent array 
     155      zhm_i_ini(1) = hginn ; zhm_i_ini(2) = hgins  ! ice thickness 
     156      zat_i_ini(1) = aginn ; zat_i_ini(2) = agins  ! ice concentration 
     157      zvt_i_ini(:) = zhm_i_ini(:) * zat_i_ini(:)   ! ice volume 
     158      zhm_s_ini(1) = hninn ; zhm_s_ini(2) = hnins  ! snow depth 
     159      zsm_i_ini(1) = sinn  ; zsm_i_ini(2) = sins   ! bulk ice salinity 
     160 
     161      !--------------------------------------------------------------------- 
     162      ! 3.2) Distribute ice concentration and thickness into the categories 
     163      !--------------------------------------------------------------------- 
     164      ! a gaussian distribution for ice concentration is used 
     165      ! then we check whether the distribution fullfills 
     166      ! volume and area conservation, positivity and ice categories bounds 
     167      DO i_hemis = 1, 2 
     168 
     169      ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 
     170 
     171      ! note for the great nemo engineers:  
     172      ! only very few of the WRITE statements are necessary for the reference version 
     173      ! they were one day useful, but now i personally doubt of their 
     174      ! potential for bringing anything useful 
     175 
     176      DO i_fill = jpl, 1, -1 
     177         IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 
     178            !---------------------------- 
     179            ! fill the i_fill categories 
     180            !---------------------------- 
     181            ! *** 1 category to fill 
     182            IF ( i_fill .EQ. 1 ) THEN 
     183               zht_i_ini(1,i_hemis)       = zhm_i_ini(i_hemis) 
     184               za_i_ini(1,i_hemis)        = zat_i_ini(i_hemis) 
     185               zht_i_ini(2:jpl,i_hemis)   = 0._wp 
     186               za_i_ini(2:jpl,i_hemis)    = 0._wp 
     187            ELSE 
     188 
     189            ! *** >1 categores to fill 
     190            !--- Ice thicknesses in the i_fill - 1 first categories 
     191               DO jl = 1, i_fill - 1 
     192                  zht_i_ini(jl,i_hemis)    = 0.5 * ( hi_max(jl) + hi_max(jl-1) ) 
     193               END DO 
     194 
     195            !--- jl0: most likely index where cc will be maximum 
     196               DO jl = 1, jpl 
     197                  IF ( ( zhm_i_ini(i_hemis) .GT. hi_max(jl-1) ) .AND. & 
     198                       ( zhm_i_ini(i_hemis) .LE. hi_max(jl)   ) ) THEN 
     199                     jl0 = jl 
     200                  ENDIF 
     201               END DO 
     202               jl0 = MIN(jl0, i_fill) 
     203 
     204            !--- Concentrations 
     205               za_i_ini(jl0,i_hemis)      = zat_i_ini(i_hemis) / SQRT(REAL(jpl)) 
     206               DO jl = 1, i_fill - 1 
     207                  IF ( jl .NE. jl0 ) THEN 
     208                     zsigma               = 0.5 * zhm_i_ini(i_hemis) 
     209                     zarg                 = ( zht_i_ini(jl,i_hemis) - zhm_i_ini(i_hemis) ) / zsigma 
     210                     za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2) 
     211                  ENDIF 
     212               END DO  
     213 
     214               zA = 0. ! sum of the areas in the jpl categories  
     215               DO jl = 1, i_fill - 1 
     216                 zA = zA + za_i_ini(jl,i_hemis) 
     217               END DO 
     218               za_i_ini(i_fill,i_hemis)   = zat_i_ini(i_hemis) - zA ! ice conc in the last category 
     219               IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
     220          
     221            !--- Ice thickness in the last category 
     222               zV = 0. ! sum of the volumes of the N-1 categories 
     223               DO jl = 1, i_fill - 1 
     224                  zV = zV + za_i_ini(jl,i_hemis)*zht_i_ini(jl,i_hemis) 
     225               END DO 
     226               zht_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis)  
     227               IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
     228 
     229            !--- volumes 
     230               zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zht_i_ini(:,i_hemis) 
     231               IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
     232 
     233            ENDIF ! i_fill 
     234 
     235            !--------------------- 
     236            ! Compatibility tests 
     237            !--------------------- 
     238            ! Test 1: area conservation 
     239            zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons ) 
     240            IF ( zconv .LT. 1.0e-6 ) THEN 
     241               ztest_1 = 1 
     242            ELSE  
     243              ! this write is useful 
     244              IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis)  
     245               ztest_1 = 0 
     246            ENDIF 
     247 
     248            ! Test 2: volume conservation 
     249            zV_cons = SUM(zv_i_ini(:,i_hemis)) 
     250            zconv = ABS(zvt_i_ini(i_hemis) - zV_cons) 
     251 
     252            IF ( zconv .LT. 1.0e-6 ) THEN 
     253               ztest_2 = 1 
     254            ELSE 
     255              ! this write is useful 
     256              IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 
     257                            ' zvt_i_ini = ', zvt_i_ini(i_hemis) 
     258               ztest_2 = 0 
     259            ENDIF 
     260 
     261            ! Test 3: thickness of the last category is in-bounds ? 
     262            IF ( zht_i_ini(i_fill, i_hemis) .GT. hi_max(i_fill-1) ) THEN 
     263               ztest_3 = 1 
     264            ELSE 
     265               ! this write is useful 
     266               IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zht_i_ini(i_fill,i_hemis) = ', & 
     267               zht_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 
     268               ztest_3 = 0 
     269            ENDIF 
     270 
     271            ! Test 4: positivity of ice concentrations 
     272            ztest_4 = 1 
     273            DO jl = 1, jpl 
     274               IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN  
     275                  ! this write is useful 
     276                  IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis) 
     277                  ztest_4 = 0 
     278               ENDIF 
     279            END DO 
     280 
     281         ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 
     282  
     283         ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 
     284 
     285      END DO ! i_fill 
     286 
     287      IF(lwp) THEN  
     288         WRITE(numout,*), ' ztests : ', ztests 
     289         IF ( ztests .NE. 4 ) THEN 
     290            WRITE(numout,*) 
     291            WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
     292            WRITE(numout,*), ' !!!! RED ALERT                  !!! ' 
     293            WRITE(numout,*), ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!' 
     294            WRITE(numout,*), ' !!!! Something is wrong in the LIM3 initialization procedure ' 
     295            WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
     296            WRITE(numout,*) 
     297            WRITE(numout,*), ' *** ztests is not equal to 4 ' 
     298            WRITE(numout,*), ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
     299            WRITE(numout,*), ' zat_i_ini : ', zat_i_ini(i_hemis) 
     300            WRITE(numout,*), ' zhm_i_ini : ', zhm_i_ini(i_hemis) 
     301         ENDIF ! ztests .NE. 4 
     302      ENDIF 
     303       
     304      END DO ! i_hemis 
     305 
     306      !--------------------------------------------------------------------- 
     307      ! 3.3) Space-dependent arrays for ice state variables 
     308      !--------------------------------------------------------------------- 
     309 
     310      ! Ice concentration, thickness and volume, snow depth, ice 
     311      ! salinity, ice age, surface temperature 
     312      DO jl = 1, jpl ! loop over categories 
     313         DO jj = 1, jpj 
     314            DO ji = 1, jpi 
     315               a_i(ji,jj,jl)   = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration 
     316               ht_i(ji,jj,jl)  = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj))  ! ice thickness 
     317               ht_s(ji,jj,jl)  = zidto(ji,jj) * zhm_s_ini(zhemis(ji,jj))     ! snow depth 
     318               sm_i(ji,jj,jl)  = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity 
     319               o_i(ji,jj,jl)   = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age 
     320               t_su(ji,jj,jl)  = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * t_bo(ji,jj) ! surf temp 
     321  
     322               ! ice volume, snow volume, salt content, age content 
     323               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume 
     324               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume 
     325               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 
     326               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content 
     327            END DO ! ji 
     328         END DO ! jj 
     329      END DO ! jl 
     330 
     331      ! Snow temperature and heat content 
     332      DO jk = 1, nlay_s 
     333         DO jl = 1, jpl ! loop over categories 
     334            DO jj = 1, jpj 
     335               DO ji = 1, jpi 
     336                   t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * rtt 
     337                   ! Snow energy of melting 
     338                   e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
     339                   ! Change dimensions 
     340                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
     341                   ! Multiply by volume, so that heat content in 10^9 Joules 
     342                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 
     343               END DO ! ji 
     344            END DO ! jj 
     345         END DO ! jl 
     346      END DO ! jk 
     347 
     348      ! Ice salinity, temperature and heat content 
     349      DO jk = 1, nlay_i 
     350         DO jl = 1, jpl ! loop over categories 
     351            DO jj = 1, jpj 
     352               DO ji = 1, jpi 
     353                   t_i(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1._wp - zidto(ji,jj) ) * rtt  
     354                   s_i(ji,jj,jk,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min 
     355                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
     356 
     357                   ! heat content per unit volume 
     358                   e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
     359                      +   lfus    * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 
     360                      -   rcp     * ( ztmelts - rtt ) ) 
     361 
     362                   ! Correct dimensions to avoid big values 
     363                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
     364 
     365                   ! Mutliply by ice volume, and divide by number of layers  
     366                   ! to get heat content in 10^9 J 
     367                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i 
     368               END DO ! ji 
     369            END DO ! jj 
     370         END DO ! jl 
     371      END DO ! jk 
     372 
     373      !-------------------------------------------------------------------- 
     374      ! 4) Global ice variables for output diagnostics                    |  
     375      !-------------------------------------------------------------------- 
     376      fsbbq (:,:)     = 0._wp 
     377      u_ice (:,:)     = 0._wp 
     378      v_ice (:,:)     = 0._wp 
     379      stress1_i(:,:)  = 0._wp 
     380      stress2_i(:,:)  = 0._wp 
     381      stress12_i(:,:) = 0._wp 
     382 
     383# if defined key_coupled 
     384      albege(:,:)   = 0.8 * tms(:,:) 
     385# endif 
     386 
     387      !-------------------------------------------------------------------- 
     388      ! 5) Moments for advection 
     389      !-------------------------------------------------------------------- 
     390 
     391      sxopw (:,:) = 0._wp  
     392      syopw (:,:) = 0._wp  
     393      sxxopw(:,:) = 0._wp  
     394      syyopw(:,:) = 0._wp  
     395      sxyopw(:,:) = 0._wp 
     396 
     397      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp 
     398      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp 
     399      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp 
     400      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp 
     401      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp 
     402 
     403      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp    
     404      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp    
     405      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp    
     406      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp    
     407      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp    
     408 
     409      sxsal  (:,:,:)  = 0._wp 
     410      sysal  (:,:,:)  = 0._wp 
     411      sxxsal (:,:,:)  = 0._wp 
     412      syysal (:,:,:)  = 0._wp 
     413      sxysal (:,:,:)  = 0._wp 
     414 
     415      !-------------------------------------------------------------------- 
     416      ! 6) Lateral boundary conditions                                    |  
    492417      !-------------------------------------------------------------------- 
    493418 
    494419      DO jl = 1, jpl 
     420 
    495421         CALL lbc_lnk( a_i(:,:,jl)  , 'T', 1. ) 
    496422         CALL lbc_lnk( v_i(:,:,jl)  , 'T', 1. ) 
    … …  
    498424         CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. ) 
    499425         CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. ) 
    500          ! 
     426 
    501427         CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. ) 
    502428         CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. ) 
    … …  
    515441         a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl) 
    516442      END DO 
     443       
     444      at_i (:,:) = 0.0_wp 
     445      DO jl = 1, jpl 
     446         at_i (:,:) = at_i (:,:) + a_i (:,:,jl) 
     447      END DO 
    517448 
    518449      CALL lbc_lnk( at_i , 'T', 1. ) 
    … …  
    521452      CALL lbc_lnk( fsbbq  , 'T', 1. ) 
    522453      ! 
    523       CALL wrk_dealloc( jpm, zgfactorn, zgfactors, zhin, zhis ) 
     454      !-------------------------------------------------------------------- 
     455      ! 6) ????                                                           |  
     456      !-------------------------------------------------------------------- 
     457      tn_ice (:,:,:) = t_su (:,:,:) 
     458 
    524459      CALL wrk_dealloc( jpi, jpj, zidto ) 
    525       ! 
     460      CALL wrk_dealloc( jpi, jpj, zhemis ) 
     461      CALL wrk_dealloc( jpl,   2, zht_i_ini,  za_i_ini,  zv_i_ini ) 
     462      CALL wrk_dealloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 
     463 
    526464   END SUBROUTINE lim_istate 
    527  
    528465 
    529466   SUBROUTINE lim_istate_init 
    … …  
    533470      !! ** Purpose : Definition of initial state of the ice  
    534471      !! 
    535       !! ** Method :   Read the namiceini namelist and check the parameter  
    536       !!             values called at the first timestep (nit000) 
    537       !! 
    538       !! ** input  :   namelist namiceini 
     472      !! ** Method : Read the namiceini namelist and check the parameter  
     473      !!       values called at the first timestep (nit000) 
     474      !! 
     475      !! ** input :  
     476      !!        Namelist namiceini 
     477      !! 
     478      !! history : 
     479      !!  8.5  ! 03-08 (C. Ethe) original code  
     480      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization 
    539481      !!----------------------------------------------------------------------------- 
    540       NAMELIST/namiceini/ ttest, hninn, hginn_u, aginn_u, hginn_d, aginn_d, hnins,   & 
    541          &                hgins_u, agins_u, hgins_d, agins_d, sinn, sins 
     482      NAMELIST/namiceini/ ttest, hninn, hnins, hginn, hgins, aginn, agins, sinn, sins 
    542483      !!----------------------------------------------------------------------------- 
    543       ! 
    544       REWIND ( numnam_ice )               ! Read Namelist namiceini  
     484 
     485      ! Define the initial parameters 
     486      ! ------------------------- 
     487 
     488      ! Read Namelist namiceini  
     489      REWIND ( numnam_ice ) 
    545490      READ   ( numnam_ice , namiceini ) 
    546       ! 
    547       IF(lwp) THEN                        ! control print 
     491      IF(lwp) THEN 
    548492         WRITE(numout,*) 
    549493         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 
    … …  
    551495         WRITE(numout,*) '   threshold water temp. for initial sea-ice    ttest      = ', ttest 
    552496         WRITE(numout,*) '   initial snow thickness in the north          hninn      = ', hninn 
    553          WRITE(numout,*) '   initial undef ice thickness in the north     hginn_u    = ', hginn_u 
    554          WRITE(numout,*) '   initial undef ice concentr. in the north     aginn_u    = ', aginn_u           
    555          WRITE(numout,*) '   initial  def  ice thickness in the north     hginn_d    = ', hginn_d 
    556          WRITE(numout,*) '   initial  def  ice concentr. in the north     aginn_d    = ', aginn_d           
    557497         WRITE(numout,*) '   initial snow thickness in the south          hnins      = ', hnins  
    558          WRITE(numout,*) '   initial undef ice thickness in the north     hgins_u    = ', hgins_u 
    559          WRITE(numout,*) '   initial undef ice concentr. in the north     agins_u    = ', agins_u           
    560          WRITE(numout,*) '   initial  def  ice thickness in the north     hgins_d    = ', hgins_d 
    561          WRITE(numout,*) '   initial  def  ice concentr. in the north     agins_d    = ', agins_d           
    562          WRITE(numout,*) '   initial  ice salinity       in the north     sinn       = ', sinn 
    563          WRITE(numout,*) '   initial  ice salinity       in the south     sins       = ', sins 
     498         WRITE(numout,*) '   initial ice thickness  in the north          hginn      = ', hginn 
     499         WRITE(numout,*) '   initial ice thickness  in the south          hgins      = ', hgins 
     500         WRITE(numout,*) '   initial ice concentr.  in the north          aginn      = ', aginn 
     501         WRITE(numout,*) '   initial ice concentr.  in the north          agins      = ', agins 
     502         WRITE(numout,*) '   initial  ice salinity  in the north          sinn       = ', sinn 
     503         WRITE(numout,*) '   initial  ice salinity  in the south          sins       = ', sins 
    564504      ENDIF 
    565       ! 
     505 
    566506   END SUBROUTINE lim_istate_init 
    567507 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r4045 r4099  
    416416 
    417417               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 )   
    418                deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 
    419 !!gm faster to replace the line above with simply: 
    420 !!                deltat(ji,jj) = MAX( delta, creepl ) 
    421 !!gm end  
    422  
     418               ! MV rewriting 
     419               ! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 
     420               !!gm faster to replace the line above with simply: 
     421               !!                deltat(ji,jj) = MAX( delta, creepl ) 
     422               !!gm end   
     423               deltat(ji,jj) = delta + creepl 
     424               ! END MV 
    423425               !-Calculate stress tensor components zs1 and zs2  
    424426               !-at centre of grid cells (see section 3.5 of CICE user's guide). 
    … …  
    741743                  &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj) 
    742744 
    743                deltat(ji,jj) = SQRT(    zdd(ji,jj)*zdd(ji,jj)   &  
    744                   &                 + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 &  
    745                   &                          ) + creepl 
    746  
     745!              deltat(ji,jj) = SQRT(    zdd(ji,jj)*zdd(ji,jj)   &  
     746!                  &                 + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 &  
     747!                  &                          ) + creepl 
     748               ! MV rewriting 
     749               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 )   
     750               deltat(ji,jj) = delta + creepl 
     751               ! END MV 
     752             
    747753            ENDIF ! zdummy 
    748754 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90

    r4045 r4099  
    371371         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    372372         t_su(:,:,jl) = z2d(:,:) 
     373         tn_ice (:,:,:) = t_su (:,:,:) 
    373374      END DO 
    374375 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r4045 r4099  
    102102      INTEGER, INTENT(in) ::   kt    ! number of iteration 
    103103      ! 
    104       INTEGER  ::   ji, jj           ! dummy loop indices 
     104      INTEGER  ::   ji, jj, jl           ! dummy loop indices 
    105105      INTEGER  ::   ierr, ifvt, i1mfr, idfr           ! local integer 
    106106      INTEGER  ::   iflt, ial , iadv , ifral, ifrdv   !   -      - 
    … …  
    145145 
    146146            !   computation the solar flux at ocean surface 
    147             zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + & 
    148                  &    ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 
     147            IF (lk_cpl) THEN ! be carfeful: not being tested yet 
     148               ! original line 
     149               !zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) 
     150               ! new line to include solar penetration (not tested) 
     151               zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 
     152               DO jl = 1, jpl 
     153                  zfcm1 = zfcm1 - qsr_ice(ji,jj,jl) * a_i(ji,jj,jl) 
     154               END DO 
     155            ELSE 
     156               zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + & 
     157                    &    ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 
     158            ENDIF 
    149159            ! fstric     Solar flux transmitted trough the ice 
    150160            ! qsr        Net short wave heat flux on free ocean 
    … …  
    227237 
    228238            !  computing freshwater exchanges at the ice/ocean interface 
    229             zemp =   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                &   - fmmec(ji,jj)                                         ! snow falling when ridging 
     239            IF (lk_cpl) THEN  
     240               zemp = - emp_tot(ji,jj) + emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
     241                  &   - rdm_snw(ji,jj) / rdt_ice 
     242            ELSE 
     243               zemp =   emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction 
     244                  &   - tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean 
     245                  &   + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice 
     246                  &   - fmmec(ji,jj)                                         ! snow falling when ridging 
     247            ENDIF 
    233248 
    234249            ! mass flux at the ocean/ice interface (sea ice fraction) 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r4045 r4099  
    8484      REAL(wp) ::   zdhnm, zhnnew, zhisn, zihic, zzc       ! 
    8585      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    86       REAL(wp) ::   zds          ! increment of bottom ice salinity 
    8786      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
    8887      REAL(wp) ::   zsm_snowice  ! snow-ice salinity 
    … …  
    429428                  zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
    430429                     &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  
    431                   zds         = zfracs * sss_m(ii,ij) - s_i_new(ji) 
     430                  zfracs = MIN( 0.5 , zfracs ) 
    432431                  s_i_new(ji) = zfracs * sss_m(ii,ij) 
    433432               ENDIF ! fc_bo_i 
    … …  
    438437         DO ji = kideb, kiut 
    439438            IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN 
    440                ! New ice salinity must not exceed 15 psu 
     439               ! New ice salinity must not exceed 20 psu 
    441440               s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 
    442441               ! Metling point in K 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r4072 r4099  
    2626   USE limthd 
    2727   USE limsbc 
    28    USE limdia 
    2928   USE limdiahsb 
    3029   USE limwri 
    … …  
    402401      ! 2.11) Ice salinity 
    403402      !--------------------- 
    404       !clem@bug: smv_i should be updated too: smv_i(:,:,:) = smv_i(:,:,:) + sm_i(:,:,:) * ( v_i(:,:,:) - zviold(:,:,:) ) 
     403      ! clem correct bug on smv_i 
     404      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
     405 
    405406      IF (  num_sal == 2  ) THEN ! general case 
    406  
    407407         DO jl = 1, jpl 
    408408            !DO jk = 1, nlay_i 
    … …  
    410410                  DO ji = 1, jpi 
    411411                     ! salinity stays in bounds 
    412                      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) ) 
    413                      i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl) + epsi20)) 
    414                      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) 
     412                     !clem 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) ) 
     413                     smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 
     414                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20 ) ) 
     415                     smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 
    415416                  END DO ! ji 
    416417               END DO ! jj 
    417418            !END DO !jk 
    418419         END DO !jl 
    419  
    420420      ENDIF 
    421421 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r4072 r4099  
    2626   USE limthd 
    2727   USE limsbc 
    28    USE limdia 
    2928   USE limdiahsb 
    3029   USE limwri 
    … …  
    305304                     zesum = zesum + e_i(ji,jj,jk,jl) 
    306305                  END DO 
    307                   IF (ind_im .LT. nlay_i ) THEN 
    308                      smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) ) 
    309                   ENDIF 
     306                  !clem IF (ind_im .LT. nlay_i ) THEN 
     307                  !clem   smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) ) 
     308                  !clem ENDIF 
    310309                  ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) 
    311310                  v_i(ji,jj,jl)  = ht_i(ji,jj,jl) * a_i(ji,jj,jl) 
    … …  
    532531      ! 2.11) Ice salinity 
    533532      !--------------------- 
    534       !clem@bug: smv_i should be updated too: smv_i(:,:,:) = smv_i(:,:,:) + sm_i(:,:,:) * ( v_i(:,:,:) - zviold(:,:,:) ) 
     533      ! clem correct bug on smv_i 
     534      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
     535 
    535536      IF (  num_sal == 2  ) THEN ! general case 
    536537         DO jl = 1, jpl 
    … …  
    539540                  DO ji = 1, jpi 
    540541                     ! salinity stays in bounds 
    541                      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) ) 
    542                      i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl) + epsi20)) 
    543                      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) 
     542                     !clem 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) ) 
     543                     smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 
     544                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20 ) ) 
     545                     smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 
    544546                  END DO ! ji 
    545547               END DO ! jj 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r3625 r4099  
    3030 
    3131   PUBLIC sbc_ice_alloc ! called in iceini(_2).F90 
     32 
     33   CHARACTER (len=8), PUBLIC :: cn_iceflx = 'none'                !: Flux handling over ice categories 
     34   LOGICAL, PUBLIC :: ln_iceflx_ave    = .FALSE. ! Average heat fluxes over all ice categories 
     35   LOGICAL, PUBLIC :: ln_iceflx_linear = .FALSE. ! Redistribute mean heat fluxes over all ice categories, using ice temperature and albedo 
    3236 
    3337# if defined  key_lim2 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r3914 r4099  
    471471      ! Coupled case: since cloud cover is not received from atmosphere  
    472472      !               ===> defined as constant value -> definition done in sbc_cpl_init 
    473       fr1_i0(:,:) = 0.18 
    474       fr2_i0(:,:) = 0.82 
     473      IF ( ALLOCATED (fr1_i0)) fr1_i0 (:,:) = 0.18 
     474      IF ( ALLOCATED (fr2_i0)) fr2_i0 (:,:) = 0.82 
    475475      !                                                      ! ------------------------- ! 
    476476      !                                                      !      10m wind module      !    
    … …  
    931931      CALL wrk_alloc( jpi,jpj, ztx, zty ) 
    932932 
    933       IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1    
     933!AC Pour eviter un stress nul sur la glace dans le cas mixed oce-ice 
     934      IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN   ;   itx =  jpr_itx1    
    934935      ELSE                                ;   itx =  jpr_otx1 
    935936      ENDIF 
    … …  
    938939      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN 
    939940 
    940          !                                                      ! ======================= ! 
    941          IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   ! 
    942             !                                                   ! ======================= ! 
     941         !                                                                                              ! ======================= ! 
     942!AC Pour eviter un stress nul sur la glace dans le cas mixes oce-ice 
     943         IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN               !   ice stress received   ! 
     944            !                                                                                           ! ======================= ! 
    943945            !   
    944946            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90

    r3625 r4099  
    9999          
    100100         fr_i(:,:) = tfreez( sss_m ) * tmask(:,:,1)      ! sea surface freezing temperature [Celcius] 
    101 #if defined key_coupled  
     101 
     102! OM : probleme. a_i pas defini dans les cas lim3 et cice 
     103#if defined key_coupled && defined key_lim2 
    102104         a_i(:,:,1) = fr_i(:,:)          
    103105#endif 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r4042 r4099  
    3232   USE sbcblk_core     ! Surface boundary condition: CORE bulk 
    3333   USE sbcblk_clio     ! Surface boundary condition: CLIO bulk 
     34   USE sbccpl          ! Surface boundary condition: coupled interface 
    3435   USE albedo          ! ocean & ice albedo 
    3536 
    … …  
    4243   USE limitd_me       ! Mechanics on ice thickness distribution 
    4344   USE limsbc          ! sea surface boundary condition 
    44    USE limdia          ! Ice diagnostics 
    4545   USE limdiahsb       ! Ice budget diagnostics 
    4646   USE limwri          ! Ice outputs 
    … …  
    7777   !!---------------------------------------------------------------------- 
    7878CONTAINS 
     79 
     80   FUNCTION fice_cell_ave ( ptab) 
     81      !!-------------------------------------------------------------------------- 
     82      !! * Compute average over categories, for grid cell (ice covered and free ocean) 
     83      !!-------------------------------------------------------------------------- 
     84      REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave 
     85      REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab 
     86      INTEGER :: jl ! Dummy loop index 
     87       
     88      fice_cell_ave (:,:) = 0.0_wp 
     89       
     90      DO jl = 1, jpl 
     91         fice_cell_ave (:,:) = fice_cell_ave (:,:) & 
     92            &                  + a_i (:,:,jl) * ptab (:,:,jl) 
     93      END DO 
     94       
     95   END FUNCTION fice_cell_ave 
     96    
     97   FUNCTION fice_ice_ave ( ptab) 
     98      !!-------------------------------------------------------------------------- 
     99      !! * Compute average over categories, for ice covered part of grid cell 
     100      !!-------------------------------------------------------------------------- 
     101      REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave 
     102      REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab 
     103 
     104      fice_ice_ave (:,:) = 0.0_wp 
     105      WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:) 
     106 
     107   END FUNCTION fice_ice_ave 
     108 
     109   !!====================================================================== 
    79110 
    80111   SUBROUTINE sbc_ice_lim( kt, kblk ) 
    … …  
    104135      REAL(wp) ::   zcoef   ! local scalar 
    105136      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky 
     137      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice (for coupled) 
     138 
     139      REAL(wp), POINTER, DIMENSION(:,:) :: zalb_ice_all    ! Mean albedo over all categories 
     140      REAL(wp), POINTER, DIMENSION(:,:) :: ztem_ice_all    ! Mean temperature over all categories 
     141       
     142      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_ice_all   ! Mean solar heat flux over all categories 
     143      REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_ice_all   ! Mean non solar heat flux over all categories 
     144      REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_ice_all   ! Mean latent heat flux over all categories 
     145      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories 
     146      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all  ! Mean d(qla)/dT over all categories 
    106147      !!---------------------------------------------------------------------- 
    107148 
     149      !- O.M. : why do we allocate all these arrays even when MOD( kt-1, nn_fsbc ) /= 0 ????? 
     150 
    108151      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
    109152 
    110153      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
     154 
     155      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
     156         CALL wrk_alloc( jpi,jpj,jpl, zalb_ice) 
     157      END IF 
     158      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
     159         CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
     160      ENDIF 
     161 
    111162 
    112163      IF( kt == nit000 ) THEN 
    … …  
    139190            t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) ) 
    140191         END DO 
     192 
     193         IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) ) 
     194          
     195         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
     196            ! 
     197            ! Compute mean albedo and temperature 
     198            zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) )  
     199            ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) )  
     200            ! 
     201         ENDIF 
    141202                                                     ! Bulk formulea - provides the following fields: 
    142203         ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2] 
    … …  
    161222               &                      tprecip   , sprecip   ,                            & 
    162223               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  ) 
     224            ! 
     225         CASE ( 5 ) 
     226            zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) ) 
     227             
     228            CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) 
     229 
     230            CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=tn_ice    ) 
     231 
     232            ! Latent heat flux is forced to 0 in coupled : 
     233            !  it is included in qns (non-solar heat flux) 
     234            qla_ice  (:,:,:) = 0.0e0_wp 
     235            dqla_ice (:,:,:) = 0.0e0_wp 
     236            ! 
    163237         END SELECT 
     238 
     239         ! Average over all categories 
     240         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
     241 
     242            z_qns_ice_all  (:,:) = fice_ice_ave ( qns_ice  (:,:,:) ) 
     243            z_qsr_ice_all  (:,:) = fice_ice_ave ( qsr_ice  (:,:,:) ) 
     244            z_dqns_ice_all (:,:) = fice_ice_ave ( dqns_ice (:,:,:) ) 
     245            z_qla_ice_all  (:,:) = fice_ice_ave ( qla_ice  (:,:,:) ) 
     246            z_dqla_ice_all (:,:) = fice_ice_ave ( dqla_ice (:,:,:) ) 
     247 
     248            DO jl = 1, jpl 
     249               dqns_ice (:,:,jl) = z_dqns_ice_all (:,:) 
     250               dqla_ice (:,:,jl) = z_dqla_ice_all (:,:) 
     251            END DO 
     252            ! 
     253            IF ( ln_iceflx_ave ) THEN 
     254               DO jl = 1, jpl 
     255                  qns_ice  (:,:,jl) = z_qns_ice_all  (:,:) 
     256                  qsr_ice  (:,:,jl) = z_qsr_ice_all  (:,:) 
     257                  qla_ice  (:,:,jl) = z_qla_ice_all  (:,:) 
     258               END DO 
     259            END IF 
     260            ! 
     261            IF ( ln_iceflx_linear ) THEN 
     262               DO jl = 1, jpl 
     263                  qns_ice  (:,:,jl) = z_qns_ice_all(:,:) + z_dqns_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:)) 
     264                  qla_ice  (:,:,jl) = z_qla_ice_all(:,:) + z_dqla_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:)) 
     265                  qsr_ice  (:,:,jl) = (1.0e0_wp-zalb_ice(:,:,jl)) / (1.0e0_wp-zalb_ice_all(:,:)) * z_qsr_ice_all(:,:) 
     266               END DO 
     267            END IF 
     268         END IF 
    164269 
    165270         !                                           !----------------------! 
    … …  
    264369         ! 
    265370         !                                           ! Diagnostics and outputs  
    266          IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   & 
    267             &             CALL lim_dia  
    268          IF (ln_limdiahsb) CALL lim_diahsb 
     371         IF (ln_limdiaout) CALL lim_diahsb 
     372!clem # if ! defined key_iomput 
    269373                          CALL lim_wri( 1  )              ! Ice outputs  
     374!clem # endif 
    270375         IF( kt == nit000 )   CALL iom_close( numrir )  ! clem: close input ice restart file 
    271376         IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file  
    … …  
    287392      ! 
    288393      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
     394      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
     395         CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice) 
     396      END IF 
     397      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
     398         CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
     399      ENDIF 
    289400      ! 
    290401      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim') 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4054 r4099  
    8484      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   & 
    8585         &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   & 
    86          &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw 
     86         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, cn_iceflx 
    8787      !!---------------------------------------------------------------------- 
    8888 
    … …  
    117117         WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs 
    118118         WRITE(numout,*) '              coupled    formulation (T if key_sbc_cpl)  ln_cpl      = ', ln_cpl 
     119         WRITE(numout,*) '              Flux handling over ice categories          cn_iceflx   = ', TRIM (cn_iceflx) 
    119120         WRITE(numout,*) '           Misc. options of sbc : ' 
    120121         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn 
    … …  
    128129      ENDIF 
    129130 
     131      !   Flux handling over ice categories  
     132      SELECT CASE ( TRIM (cn_iceflx)) 
     133      CASE ('ave') 
     134         ln_iceflx_ave    = .TRUE. 
     135         ln_iceflx_linear = .FALSE. 
     136      CASE ('linear') 
     137         ln_iceflx_ave    = .FALSE. 
     138         ln_iceflx_linear = .TRUE. 
     139      CASE default 
     140         ln_iceflx_ave    = .FALSE. 
     141         ln_iceflx_linear = .FALSE. 
     142      END SELECT 
     143      IF(lwp) WRITE(numout,*) '              Fluxes averaged over all ice categories         ln_iceflx_ave    = ', ln_iceflx_ave 
     144      IF(lwp) WRITE(numout,*) '              Fluxes distributed linearly over ice categories ln_iceflx_linear = ', ln_iceflx_linear 
     145      ! 
    130146      !                              ! allocate sbc arrays 
    131147      IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' ) 
    … …  
    166182      IF( ( nn_ice == 3 .OR. nn_ice == 4 ) .AND. nn_ice_embd == 0 )   & 
    167183         &   CALL ctl_stop( 'LIM3 and CICE sea-ice models require nn_ice_embd = 1 or 2' ) 
     184 
     185      IF( ln_iceflx_ave .AND. ln_iceflx_linear ) & 
     186         &   CALL ctl_stop( ' ln_iceflx_ave and ln_iceflx_linear options are not compatible' ) 
     187 
     188      IF( ( nn_ice ==3 .AND. lk_cpl) .AND. .NOT. ( ln_iceflx_ave .OR. ln_iceflx_linear ) ) & 
     189         &   CALL ctl_stop( ' With lim3 coupled, either ln_iceflx_ave or ln_iceflx_linear must be set to .TRUE.' ) 
    168190       
    169191      IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag 
    … …  
    295317      CASE(  2 )   ;         CALL sbc_ice_lim_2( kt, nsbc )          ! LIM-2 ice model 
    296318      CASE(  3 )   ;         CALL sbc_ice_lim  ( kt, nsbc )          ! LIM-3 ice model 
     319      !is it useful? 
     320      !CASE(  4 )   ;         CALL sbc_ice_cice ( kt, nsbc )          ! CICE ice model 
    297321      END SELECT                                               
    298322 
Note: See TracChangeset for help on using the changeset viewer.