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 4332 for branches – NEMO

Changeset 4332 for branches


Ignore:
Timestamp:
2013-12-11T15:38:42+01:00 (10 years ago)
Author:
clem
Message:

update LIM3 to fix remaining bugs. Now working in global and regional config.

Location:
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM
Files:
1 deleted
30 edited

Legend:

Unmodified
Added
Removed

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

    r4099 r4332  
    77ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    88AMM12 OPA_SRC 
     9ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 
    910ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    10 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 

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

    r4220 r4332  
    292292   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dh_i_surf2D, dh_i_bott2D, fstbif, fsup2D, focea2D, q_s 
    293293 
    294    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   patho_case ! number of the pathological case (if any, of course) 
    295  
    296294   !!-------------------------------------------------------------------------- 
    297295   !! * Ice global state variables 
    … …  
    396394   LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F) 
    397395   LOGICAL               , PUBLIC ::   ln_nicep      = .FALSE.             !: flag for sea-ice points output (T) or not (F) 
    398    REAL(wp)              , PUBLIC ::   cai           = 1.40e-3            !: atmospheric drag over sea ice 
    399    REAL(wp)              , PUBLIC ::   cao           = 1.00e-3            !: atmospheric drag over ocean 
     396   REAL(wp)              , PUBLIC ::   cai           = 1.4e-3            !: atmospheric drag over sea ice 
     397   REAL(wp)              , PUBLIC ::   cao           = 1.0e-3            !: atmospheric drag over ocean 
    400398   REAL(wp)              , PUBLIC ::   amax          = 0.99               !: maximum ice concentration 
    401399   ! 
    … …  
    465463         &      fsup2D     (jpi,jpj) , focea2D    (jpi,jpj) , q_s   (jpi,jpj) , STAT=ierr(ii) ) 
    466464 
    467       ii = ii + 1 
    468       ALLOCATE( patho_case(jpi, jpj, jpl)  , STAT=ierr(ii) ) 
    469  
    470465      ! * Ice global state variables 
    471466      ii = ii + 1 

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

    r4220 r4332  
    8888      ! 
    8989      !                                ! Initial sea-ice state 
    90       IF( .NOT.ln_rstart ) THEN              ! start from rest 
     90      IF( .NOT. ln_rstart ) THEN              ! start from rest 
    9191         numit = 0 
    9292         numit = nit000 - 1 
    … …  
    101101      ENDIF 
    102102      ! 
    103       CALL lim_sbc_init                ! ice surface boundary condition    
    104       ! 
    105       fr_i(:,:) = at_i(:,:)           ! initialisation of sea-ice fraction 
    106       tn_ice(:,:,:) = t_su(:,:,:) 
     103      hi_max(jpl) = 99._wp  ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl) 
     104      ! 
     105      CALL lim_sbc_init                 ! ice surface boundary condition    
     106      ! 
     107      fr_i(:,:) = at_i(:,:)             ! initialisation of sea-ice fraction 
     108      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu 
    107109      ! 
    108110      nstart = numit  + nn_fsbc       
    … …  
    132134      READ  ( numnam_ice , namicerun ) 
    133135      ! 
    134       IF( lk_mpp .AND. ln_nicep ) THEN 
    135          ln_nicep = .FALSE. 
    136          CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' ) 
    137       ENDIF 
     136      !IF( lk_mpp .AND. ln_nicep ) THEN 
     137      !   ln_nicep = .FALSE. 
     138      !   CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' ) 
     139      !ENDIF 
    138140      ! 
    139141      IF(lwp) THEN                        ! control print 
    … …  
    159161      !! ** Purpose :   Initializes the ice thickness distribution 
    160162      !! ** Method  :   ... 
     163      !!    Note    : hi_max(jpl) is here set up to a value close to 7 m for 
     164      !!              limistate (only) and is changed to 99 m in ice_init 
    161165      !!------------------------------------------------------------------ 
    162166      INTEGER  ::   jl, jm               ! dummy loop index 

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

    r4099 r4332  
    4444      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables 
    4545      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables 
    46       REAL(wp) :: epsi06 = 1.0e-6 
    47    
     46      REAL(wp)                     ::   epsi06 = 1.0e-6 
     47      REAL(wp)                     ::   zc1, zc2, zc3, zx1, zdh   ! local scalars 
     48      REAL(wp), DIMENSION(0:jpl)   ::   zhi_max         !:Boundary of ice thickness categories in thickness space 
     49  
    4850      IF( nn_timing == 1 )  CALL timing_start('limcat_1D') 
    4951      !-------------------------------------------------------------------- 
    … …  
    5153      !-------------------------------------------------------------------- 
    5254      ijpij = SIZE(zhti,1) 
    53       zht_i(1:ijpij,1:jpl) = 0.d0 
    54       zht_s(1:ijpij,1:jpl) = 0.d0 
    55       za_i (1:ijpij,1:jpl) = 0.d0 
     55      zht_i(1:ijpij,1:jpl) = 0._wp 
     56      zht_s(1:ijpij,1:jpl) = 0._wp 
     57      za_i (1:ijpij,1:jpl) = 0._wp 
    5658 
    5759      !------------------------------------------------------------------------------------ 
    … …  
    6668      !         fulfills ice volume concervation between input and output (ztests=4)  
    6769      !-------------------------------------------------------------------------------------- 
     70 
     71      !- Thickness categories boundaries 
     72      !    hi_max is calculated in iceini.F90 but since limcat_1D.F90 routine  
     73      !    is called before (in bdydta.F90), one must recalculate it.    
     74      !    Note clem: there may be a way of doing things cleaner 
     75      !---------------------------------- 
     76      zhi_max(:) = 0._wp 
     77      zc1 =  3._wp / REAL( jpl , wp ) ; zc2 = 10._wp * zc1 ; zc3 =  3._wp 
     78      DO jl = 1, jpl 
     79         zx1 = REAL( jl-1 , wp ) / REAL( jpl , wp ) 
     80         zhi_max(jl) = zhi_max(jl-1) + zc1 + zc2 * ( 1._wp + TANH( zc3 * ( zx1 - 1._wp ) ) ) 
     81      END DO 
    6882  
    6983      ! ---------------------------------------- 
    … …  
    7185      ! ---------------------------------------- 
    7286      DO ji = 1, ijpij 
    73         ! snow thickness in each category  
    74          zht_s(ji,1:jpl) = zhts(ji) 
    7587          
     88         IF( zhti(ji) > 0._wp ) THEN 
     89 
    7690         ! initialisation of tests 
    7791         ztest_1 = 0 
    … …  
    87101             
    88102            ! initialisation of ice variables for each try 
    89             zht_i(ji,1:jpl) = 0.d0 
    90             za_i (ji,1:jpl) = 0.d0 
     103            zht_i(ji,1:jpl) = 0._wp 
     104            za_i (ji,1:jpl) = 0._wp 
    91105             
    92106            ! *** case very thin ice: fill only category 1 
    … …  
    94108               zht_i(ji,1) = zhti(ji) 
    95109               za_i (ji,1) = zai (ji) 
    96                 
    97110               ! *** case ice is thicker: fill categories >1 
    98111            ELSE 
    … …  
    100113               ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2  
    101114               DO jl = 1, i_fill - 1 
    102                   zht_i(ji,jl) = ( hi_max(jl) + hi_max(jl-1) ) / 2. 
     115                  zht_i(ji,jl) = ( zhi_max(jl) + zhi_max(jl-1) ) * 0.5_wp 
    103116               END DO 
    104117                
    … …  
    106119               jl0 = i_fill 
    107120               DO jl = 1, i_fill 
    108                   IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
     121                  IF ( ( zhti(ji) >= zhi_max(jl-1) ) .AND. ( zhti(ji) < zhi_max(jl) ) ) THEN 
    109122                     jl0 = jl 
    110123           CYCLE 
    … …  
    116129               DO jl = 1, i_fill - 1 
    117130                  IF ( jl == jl0 ) CYCLE 
    118                   zarg           = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) / 2. ) 
     131                  zarg           = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
    119132                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
    120133               END DO 
    … …  
    141154             
    142155            ! Test 3: thickness of the last category is in-bounds ? 
    143             IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) ztest_3 = 1 
     156            IF ( zht_i(ji,i_fill) >= zhi_max(i_fill-1) ) ztest_3 = 1 
    144157             
    145158            ! Test 4: positivity of ice concentrations 
    146159            ztest_4 = 1 
    147160            DO jl = 1, i_fill 
    148                IF ( za_i(ji,jl) < 0.0d0 ) ztest_4 = 0 
     161               IF ( za_i(ji,jl) < 0._wp ) ztest_4 = 0 
    149162            END DO 
    150163             
    … …  
    164177         !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:) 
    165178         !ENDIF 
    166              
     179          
     180         ENDIF ! if zhti > 0 
     181 
    167182      END DO ! i loop 
     183 
     184      ! ------------------------------------------------ 
     185      ! Adding Snow in each category where za_i is not 0 
     186      ! ------------------------------------------------  
     187      DO jl = 1, jpl 
     188         DO ji = 1, ijpij 
     189            IF( za_i(ji,jl) > 0._wp ) THEN 
     190               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 
     191               ! In case snow load is in excess that would lead to transformation from snow to ice 
     192               ! Then, transfer the snow excess into the ice (different from limthd_dh) 
     193               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )  
     194               ! recompute ht_i, ht_s avoiding out of bounds values 
     195               zht_i(ji,jl) = MIN( zhi_max(jl), zht_i(ji,jl) + zdh ) 
     196               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn ) 
     197            ENDIF 
     198         ENDDO 
     199      ENDDO 
    168200 
    169201      IF( nn_timing == 1 )  CALL timing_stop('limcat_1D') 

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

    r4045 r4332  
    2929 
    3030   LOGICAL  ::   linit = .TRUE.              ! initialization flag (set to flase after the 1st call) 
    31    REAL(wp) ::   epsi04 = 1e-04              ! constant 
     31   REAL(wp) ::   epsi04 = 1.e-04              ! constant 
    3232   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient 
    3333 
    … …  
    5454      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied 
    5555      ! 
    56       INTEGER  ::  ji, jj            ! dummy loop indices 
    57       INTEGER  ::  its, iter, ierr   ! local integers 
    58       REAL(wp) ::   zalfa, zrlxint, zconv, zeps   ! local scalars 
     56      INTEGER  ::  ji, jj                   ! dummy loop indices 
     57      INTEGER  ::  its, iter, ierr          ! local integers 
     58      REAL(wp) ::   zalfa, zrlxint, zconv   ! local scalars 
    5959      REAL(wp), POINTER, DIMENSION(:,:) ::   zrlx, zflu, zflv, zdiv0, zdiv, ztab0 
    6060      CHARACTER(lc) ::   charout   ! local character 
    … …  
    7979      zalfa = 0.5_wp                      ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 
    8080      its   = 100                         ! Maximum number of iteration 
    81       zeps  =  2._wp * epsi04 
    8281      ! 
    8382      ztab0(:, : ) = ptab(:,:)      ! Arrays initialization 
    … …  
    9493      iter  = 0 
    9594      ! 
    96       DO WHILE( zconv > zeps .AND. iter <= its )   ! Sub-time step loop 
     95      DO WHILE( zconv > ( 2._wp * epsi04 ) .AND. iter <= its )   ! Sub-time step loop 
    9796         ! 
    9897         iter = iter + 1                                 ! incrementation of the sub-time step number 

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

    r4099 r4332  
    7676      !!                6) Lateral boundary conditions 
    7777      !! 
     78      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even 
     79      !!              where there is no ice (clem: I do not know why but it is mandatory)  
     80      !! 
    7881      !! History : 
    7982      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code 
    … …  
    8487      !! * Local variables 
    8588      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices 
    86       REAL(wp)   :: epsi06, epsi20, ztmelts 
     89      REAL(wp)   :: epsi20, ztmelts, zdh 
    8790      INTEGER    :: i_hemis, i_fill, jl0   
    8891      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 
    … …  
    98101      CALL wrk_alloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 
    99102 
    100       epsi06   = 1.0e-6 
    101103      epsi20   = 1.0e-20 
    102104      IF(lwp) WRITE(numout,*) 
    … …  
    308310      !--------------------------------------------------------------------- 
    309311 
    310       ! Ice concentration, thickness and volume, snow depth, ice 
    311       ! salinity, ice age, surface temperature 
     312      ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 
    312313      DO jl = 1, jpl ! loop over categories 
    313314         DO jj = 1, jpj 
    … …  
    315316               a_i(ji,jj,jl)   = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration 
    316317               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 
     318               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zhm_s_ini( zhemis(ji,jj) ) / zhm_i_ini( zhemis(ji,jj) ) )  ! snow depth 
     319               sm_i(ji,jj,jl)  = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min ! salinity 
    319320               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 
     321               t_su(ji,jj,jl)  = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * 270.0 ! surf temp 
     322 
     323               ! This case below should not be used if (ht_s/ht_i) is ok in namelist 
     324               ! In case snow load is in excess that would lead to transformation from snow to ice 
     325               ! Then, transfer the snow excess into the ice (different from limthd_dh) 
     326               zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 )  
     327               ! recompute ht_i, ht_s avoiding out of bounds values 
     328               ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh ) 
     329               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn ) 
     330 
     331               ! ice volume, salt content, age content 
    323332               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume 
    324333               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume 

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

    r4072 r4332  
    4343   PUBLIC   lim_itd_me_alloc        ! called by iceini.F90 
    4444 
    45    REAL(wp) ::   epsi11 = 1.e-11_wp   ! constant values 
     45   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values 
    4646   REAL(wp) ::   epsi10 = 1.e-10_wp   ! constant values 
    4747   REAL(wp) ::   epsi06 = 1.e-06_wp   ! constant values 
    … …  
    5050   ! Variables shared among ridging subroutines 
    5151   !----------------------------------------------------------------------- 
    52    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
    53    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
    54    ! 
    55    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/ 
    56    !                                                           !  closing associated w/ category n 
    57    ! 
    58    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
    59    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
    60    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice 
    61    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice  
    62    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging 
    63    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting 
     52   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area 
     53   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged 
     54   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/ 
     55   !                                                     ! closing associated w/ category n 
     56   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness 
     57   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness 
     58   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice 
     59   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice  
     60   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging 
     61   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting 
    6462 
    6563   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier 
    6664   REAL(wp), PARAMETER ::   kraft   = 2.0_wp    ! rafting multipliyer 
    67    REAL(wp), PARAMETER ::   kamax   = 1.0 
     65   REAL(wp), PARAMETER ::   kamax   = 1.0_wp    ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99) 
    6866 
    6967   REAL(wp) ::   Cp                             !  
    … …  
    183181      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    184182      !-----------------------------------------------------------------------------! 
    185       ! Set hi_max(ncat) to a big value to ensure that all ridged ice is thinner than hi_max(ncat). 
    186  
    187       hi_max(jpl) = 999.99 
    188  
    189183      Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0                ! proport const for PE 
    190184      ! 
    … …  
    265259               ! Reduce the closing rate if more than 100% of the open water  
    266260               ! would be removed.  Reduce the opening rate proportionately. 
    267                IF ( ato_i(ji,jj) .GT. epsi11 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 
     261               IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 
    268262                  w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 
    269263                  IF ( w1 .GT. ato_i(ji,jj)) THEN 
    … …  
    285279            DO jj = 1, jpj 
    286280               DO ji = 1, jpi 
    287                   IF ( a_i(ji,jj,jl) > epsi11 .AND. athorn(ji,jj,jl) > 0._wp )THEN 
     281                  IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN 
    288282                     w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
    289283                     IF ( w1  >  a_i(ji,jj,jl) ) THEN 
    … …  
    316310         DO jj = 1, jpj 
    317311            DO ji = 1, jpi 
    318                IF (ABS(asum(ji,jj) - kamax ) .LT. epsi11) THEN 
     312               IF (ABS(asum(ji,jj) - kamax ) .LT. epsi10) THEN 
    319313                  closing_net(ji,jj) = 0._wp 
    320314                  opning     (ji,jj) = 0._wp 
    … …  
    358352         DO ji = 1, jpi 
    359353 
    360             IF(ABS(asum(ji,jj) - kamax) > epsi11 ) asum_error = .true. 
     354            IF(ABS(asum(ji,jj) - kamax) > epsi10 ) asum_error = .true. 
    361355 
    362356            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice 
    … …  
    377371      DO jj = 1, jpj 
    378372         DO ji = 1, jpi 
    379             IF( ABS( asum(ji,jj) - kamax)  >  epsi11 ) THEN   ! there is a bug 
     373            IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug 
    380374               WRITE(numout,*) ' ' 
    381375               WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
    … …  
    400394 
    401395      !-----------------------------------------------------------------------------! 
    402       ! 6) Updating state variables and trend terms 
     396      ! 6) Updating state variables and trend terms (done in limupdate) 
    403397      !-----------------------------------------------------------------------------! 
    404  
    405398      CALL lim_var_glo2eqv 
    406399      CALL lim_itd_me_zapsmall 
    … …  
    419412         END DO 
    420413      END DO 
    421  
    422       !----------------- 
    423       !  Trend terms 
    424       !----------------- 
    425  
    426       d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:) 
    427       d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:) 
    428       d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:) 
    429       d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)   
    430       d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)    
    431       d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)   
    432       d_e_i_trp  (:,:,:,:) = e_i  (:,:,:,:) - old_e_i  (:,:,:,:) 
    433       d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
    434       d_smv_i_trp(:,:,:)   = 0._wp 
    435       IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
    436414 
    437415      IF(ln_ctl) THEN     ! Control print 
    … …  
    491469      ! ------------------------------- 
    492470 
    493       !-------------------------! 
    494       ! Back to initial values 
    495       !-------------------------! 
    496  
    497       ! update of fields will be made later in lim update 
    498       u_ice(:,:)   = old_u_ice(:,:) 
    499       v_ice(:,:)   = old_v_ice(:,:) 
    500       a_i(:,:,:)   = old_a_i(:,:,:) 
    501       v_s(:,:,:)   = old_v_s(:,:,:) 
    502       v_i(:,:,:)   = old_v_i(:,:,:) 
    503       e_s(:,:,:,:) = old_e_s(:,:,:,:) 
    504       e_i(:,:,:,:) = old_e_i(:,:,:,:) 
    505       oa_i(:,:,:)  = old_oa_i(:,:,:) 
    506       IF(  num_sal == 2  )   smv_i(:,:,:) = old_smv_i(:,:,:) 
    507  
    508       !----------------------------------------------------! 
    509       ! Advection of ice in a free cell, newly ridged ice 
    510       !----------------------------------------------------! 
    511  
    512       ! to allow for thermodynamics to melt new ice 
    513       ! we immediately advect ice in free cells 
    514  
    515       ! heat content has to be corrected before ice volume 
    516 !clem@order 
    517 !      DO jl = 1, jpl 
    518 !         DO jk = 1, nlay_i 
    519 !            DO jj = 1, jpj 
    520 !               DO ji = 1, jpi 
    521 !                  IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 
    522 !                     ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 
    523 !                     old_e_i(ji,jj,jk,jl)   = d_e_i_trp(ji,jj,jk,jl) 
    524 !                     d_e_i_trp(ji,jj,jk,jl) = 0._wp 
    525 !                  ENDIF 
    526 !               END DO 
    527 !            END DO 
    528 !         END DO 
    529 !      END DO 
    530 ! 
    531 !      DO jl = 1, jpl 
    532 !         DO jj = 1, jpj 
    533 !            DO ji = 1, jpi 
    534 !               IF(  old_v_i  (ji,jj,jl) < epsi06  .AND. & 
    535 !                    d_v_i_trp(ji,jj,jl) > epsi06    ) THEN 
    536 !                  old_v_i   (ji,jj,jl)   = d_v_i_trp(ji,jj,jl) 
    537 !                  d_v_i_trp (ji,jj,jl)   = 0._wp 
    538 !                  old_a_i   (ji,jj,jl)   = d_a_i_trp(ji,jj,jl) 
    539 !                  d_a_i_trp (ji,jj,jl)   = 0._wp 
    540 !                  old_v_s   (ji,jj,jl)   = d_v_s_trp(ji,jj,jl) 
    541 !                  d_v_s_trp (ji,jj,jl)   = 0._wp 
    542 !                  old_e_s   (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 
    543 !                  d_e_s_trp (ji,jj,1,jl) = 0._wp 
    544 !                  old_oa_i  (ji,jj,jl)   = d_oa_i_trp(ji,jj,jl) 
    545 !                  d_oa_i_trp(ji,jj,jl)   = 0._wp 
    546 !                  IF(  num_sal == 2  )   old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 
    547 !                  d_smv_i_trp(ji,jj,jl)  = 0._wp 
    548 !               ENDIF 
    549 !            END DO 
    550 !         END DO 
    551 !      END DO 
    552 !clem@order 
    553471      ENDIF  ! ln_limdyn=.true. 
    554472      ! 
    … …  
    605523               DO ji = 1, jpi 
    606524                  ! 
    607                   IF(  a_i(ji,jj,jl)    > epsi11  .AND.     & 
    608                        athorn(ji,jj,jl) > 0._wp   ) THEN 
     525                  IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN 
    609526                     hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    610527                     !---------------------------- 
    … …  
    624541                        * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )    
    625542!!gm Optimization:  (a**3-b**3)/(a-b) = a*a+ab+b*b   ==> less costly operations even if a**3 is replaced by a*a*a...                     
    626                   ENDIF            ! aicen > epsi11 
     543                  ENDIF            ! aicen > epsi10 
    627544                  ! 
    628545               END DO ! ji 
    … …  
    681598         DO jj = 2, jpj - 1 
    682599            DO ji = 2, jpi - 1 
    683                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is 
     600               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is 
    684601                  ! present 
    685602                  zworka(ji,jj) = 4.0 * strength(ji,jj)              & 
    … …  
    721638         DO jj = 1, jpj - 1 
    722639            DO ji = 1, jpi - 1 
    723                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN       ! ice is present 
     640               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN       ! ice is present 
    724641                  numts_rm = 1 ! number of time steps for the running mean 
    725642                  IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 
    … …  
    794711      DO jj = 1, jpj 
    795712         DO ji = 1, jpi 
    796             IF( ato_i(ji,jj) > epsi11 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj) 
     713            IF( ato_i(ji,jj) > epsi10 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj) 
    797714            ELSE                               ;   Gsum(ji,jj,0) = 0._wp 
    798715            ENDIF 
    … …  
    804721         DO jj = 1, jpj  
    805722            DO ji = 1, jpi 
    806                IF( a_i(ji,jj,jl) .GT. epsi11 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 
     723               IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 
    807724               ELSE                                   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 
    808725               ENDIF 
    … …  
    887804      IF ( raftswi == 1 ) THEN 
    888805 
    889          IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi11 ) THEN 
     806         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN 
    890807            DO jl = 1, jpl 
    891808               DO jj = 1, jpj 
    892809                  DO ji = 1, jpi 
    893                      IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. & 
    894                         epsi11 ) THEN 
     810                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN 
    895811                        WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 
    896812                        WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl 
    … …  
    938854            DO ji = 1, jpi 
    939855 
    940                IF (a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 
     856               IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 
    941857                  hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    942858                  hrmean          = MAX(SQRT(Hstar*hi), hi*krdgmin) 
    … …  
    992908      INTEGER ::   ij                ! horizontal index, combines i and j loops 
    993909      INTEGER ::   icells            ! number of cells with aicen > puny 
    994       REAL(wp) ::   zeps, zindb, zsrdg2   ! local scalar 
     910      REAL(wp) ::   zindb, zsrdg2   ! local scalar 
    995911      REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration 
    996912 
    … …  
    1044960 
    1045961      IF( con_i ) THEN 
    1046          CALL lim_column_sum (jpl,   v_i, vice_init ) 
    1047          WRITE(numout,*) ' vice_init  : ', vice_init(jiindx,jjindx) 
     962         CALL lim_column_sum        (jpl,    v_i,       vice_init ) 
    1048963         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_init ) 
    1049          WRITE(numout,*) ' eice_init  : ', eice_init(jiindx,jjindx) 
     964         DO ji = mi0(jiindx), mi1(jiindx) 
     965            DO jj = mj0(jjindx), mj1(jjindx) 
     966               WRITE(numout,*) ' vice_init  : ', vice_init(ji,jj) 
     967               WRITE(numout,*) ' eice_init  : ', eice_init(ji,jj) 
     968            END DO 
     969         END DO 
    1050970      ENDIF 
    1051  
    1052       zeps   = 1.e-20_wp 
    1053971 
    1054972      !------------------------------------------------------------------------------- 
    … …  
    1062980            ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        & 
    1063981               &                        + opning(ji,jj)                          * rdt_ice 
    1064             IF( ato_i(ji,jj) < -epsi11 ) THEN 
     982            IF( ato_i(ji,jj) < -epsi10 ) THEN 
    1065983               neg_ato_i = .TRUE. 
    1066984            ELSEIF( ato_i(ji,jj) < 0._wp ) THEN    ! roundoff error 
    … …  
    1074992         DO jj = 1, jpj  
    1075993            DO ji = 1, jpi 
    1076                IF( ato_i(ji,jj) < -epsi11 ) THEN  
     994               IF( ato_i(ji,jj) < -epsi10 ) THEN  
    1077995                  WRITE(numout,*) ''   
    1078996                  WRITE(numout,*) 'Ridging error: ato_i < 0' 
    1079997                  WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) 
    1080                ENDIF               ! ato_i < -epsi11 
     998               ENDIF               ! ato_i < -epsi10 
    1081999            END DO 
    10821000         END DO 
    … …  
    11201038         DO jj = 1, jpj 
    11211039            DO ji = 1, jpi 
    1122                IF( aicen_init(ji,jj,jl1)  >  epsi11  .AND.  athorn(ji,jj,jl1)  >  0._wp       & 
    1123                   .AND. closing_gross(ji,jj) > 0._wp ) THEN 
     1040               IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp  & 
     1041                  &   .AND. closing_gross(ji,jj) > 0._wp ) THEN 
    11241042                  icells = icells + 1 
    11251043                  indxi(icells) = ji 
    … …  
    11581076            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 
    11591077 
    1160             IF (afrac(ji,jj) > kamax + epsi11) THEN  !riging 
     1078            IF (afrac(ji,jj) > kamax + epsi10) THEN  !riging 
    11611079               large_afrac = .true. 
    11621080            ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error 
    11631081               afrac(ji,jj) = kamax 
    11641082            ENDIF 
    1165             IF (afrft(ji,jj) > kamax + epsi11) THEN !rafting 
     1083            IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting 
    11661084               large_afrft = .true. 
    11671085            ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error 
    … …  
    12221140            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 
    12231141            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 
    1224             !clem diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 
    12251142            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 
    12261143 
    … …  
    12761193 
    12771194               ! corrected sea water salinity 
    1278                zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - zeps ) ) 
    1279                zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), zeps ) 
     1195               zindb  = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) ) 
     1196               zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 ) 
    12801197 
    12811198               ztmelts          = - tmut * zdummy + rtt 
    … …  
    13121229               ji = indxi(ij) 
    13131230               jj = indxj(ij) 
    1314                IF( afrac(ji,jj) > kamax + epsi11 ) THEN  
     1231               IF( afrac(ji,jj) > kamax + epsi10 ) THEN  
    13151232                  WRITE(numout,*) '' 
    13161233                  WRITE(numout,*) ' ardg > a_i' 
    … …  
    13241241               ji = indxi(ij) 
    13251242               jj = indxj(ij) 
    1326                IF( afrft(ji,jj) > kamax + epsi11 ) THEN  
     1243               IF( afrft(ji,jj) > kamax + epsi10 ) THEN  
    13271244                  WRITE(numout,*) '' 
    13281245                  WRITE(numout,*) ' arft > a_i' 
    … …  
    14241341         fieldid = ' v_i : limitd_me ' 
    14251342         CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid)  
    1426          WRITE(numout,*) ' vice_init  : ', vice_init(jiindx,jjindx) 
    1427          WRITE(numout,*) ' vice_final : ', vice_final(jiindx,jjindx) 
    14281343 
    14291344         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_final ) 
    14301345         fieldid = ' e_i : limitd_me ' 
    14311346         CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)  
    1432          WRITE(numout,*) ' eice_init  : ', eice_init(jiindx,jjindx) 
    1433          WRITE(numout,*) ' eice_final : ', eice_final(jiindx,jjindx) 
     1347 
     1348         DO ji = mi0(jiindx), mi1(jiindx) 
     1349            DO jj = mj0(jjindx), mj1(jjindx) 
     1350               WRITE(numout,*) ' vice_init  : ', vice_init (ji,jj) 
     1351               WRITE(numout,*) ' vice_final : ', vice_final(ji,jj) 
     1352               WRITE(numout,*) ' eice_init  : ', eice_init (ji,jj) 
     1353               WRITE(numout,*) ' eice_final : ', eice_final(ji,jj) 
     1354            END DO 
     1355         END DO 
    14341356      ENDIF 
    14351357      ! 
    … …  
    15331455 
    15341456      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace 
    1535        
     1457      REAL(wp)                          ::   zmask_glo 
    15361458!!gm      REAL(wp) ::   xtmp      ! temporary variable 
    15371459      !!------------------------------------------------------------------- 
    … …  
    15451467         ! Abort model in case of negative area. 
    15461468         !----------------------------------------------------------------- 
    1547          IF( MINVAL(a_i(:,:,jl)) .LT. -epsi11 .AND. ln_nicep ) THEN 
    1548             DO jj = 1, jpj 
    1549                DO ji = 1, jpi 
    1550                   IF ( a_i(ji,jj,jl) .LT. -epsi11 ) THEN 
    1551                      WRITE (numout,*) ' ALERTE 98 '  
    1552                      WRITE (numout,*) ' Negative ice area: ji, jj, jl: ', ji, jj,jl 
    1553                      WRITE (numout,*) ' a_i    ', a_i(ji,jj,jl) 
    1554                   ENDIF 
    1555                END DO 
    1556             END DO 
    1557          ENDIF 
    1558  
    15591469         icells = 0 
    1560          zmask  = 0._wp 
     1470         zmask(:,:)  = 0._wp 
    15611471         DO jj = 1, jpj 
    15621472            DO ji = 1, jpi 
    1563                IF(  ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0._wp   ) .OR.   & 
    1564                   & ( a_i(ji,jj,jl) .GT. 0._wp   .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) .OR.   & 
    1565                   & ( v_i(ji,jj,jl)  ==  0._wp   .AND. a_i(ji,jj,jl) .GT. 0._wp   ) .OR.   & 
    1566                   & ( v_i(ji,jj,jl) .GT. 0._wp   .AND. v_i(ji,jj,jl) .LT. 1.e-12  )      )   zmask(ji,jj) = 1._wp 
    1567             END DO 
    1568          END DO 
    1569          IF( ln_nicep )   WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean ' 
     1473               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   & 
     1474                  & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   & 
     1475                  & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   & 
     1476                  & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp 
     1477            END DO 
     1478         END DO 
     1479         zmask_glo = glob_sum(zmask) 
     1480         !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean ' 
    15701481 
    15711482         !----------------------------------------------------------------- 
    … …  
    15791490!!gm                  xtmp = xtmp * unit_fac 
    15801491                  ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
    1581                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1 - zmask(ji,jj) ) 
     1492                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) 
    15821493               END DO 
    15831494            END DO 
    … …  
    16011512!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ??????? 
    16021513 
    1603                t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) ) 
     1514               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 
    16041515 
    16051516               !----------------------------------------------------------------- 
    … …  
    16151526 
    16161527               ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj) 
    1617                a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1618                v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1619                v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1620                t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1 - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 
    1621                oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1622                smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
     1528               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
     1529               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
     1530               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
     1531               t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 
     1532               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
     1533               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    16231534               ! 
    16241535            END DO 

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

    r4155 r4332  
    4545   PUBLIC   lim_itd_shiftice 
    4646 
    47    REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values 
    48    REAL(wp) ::   epsi13 = 1.e-13_wp   ! 
    4947   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
     48   REAL(wp) ::   epsi06 = 1.e-6_wp   ! 
    5049 
    5150   !!---------------------------------------------------------------------- 
    … …  
    110109      CALL lim_var_glo2eqv    ! only for info 
    111110 
    112       !---------------------------------------------------------------------------------------- 
    113       !  4) Computation of trend terms and get back to old values       
    114       !---------------------------------------------------------------------------------------- 
    115  
    116       !- Trend terms 
    117       d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:)  
    118       d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:) 
    119       d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)   
    120       d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)  
    121       d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 
    122       !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - old_oa_i (:,:,:) 
    123       d_smv_i_thd(:,:,:) = 0._wp 
    124       IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
    125  
    126       ! diag only (clem) 
    127       dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    128  
    129       IF(ln_ctl) THEN   ! Control print 
     111     IF(ln_ctl) THEN   ! Control print 
    130112         CALL prt_ctl_info(' ') 
    131113         CALL prt_ctl_info(' - Cell values : ') 
    … …  
    183165      ! ------------------------------- 
    184166      ! 
    185       !- Recover Old values 
    186       a_i(:,:,:)   = old_a_i (:,:,:) 
    187       v_s(:,:,:)   = old_v_s (:,:,:) 
    188       v_i(:,:,:)   = old_v_i (:,:,:) 
    189       e_s(:,:,:,:) = old_e_s (:,:,:,:) 
    190       e_i(:,:,:,:) = old_e_i (:,:,:,:) 
    191       !?? oa_i(:,:,:)  = old_oa_i(:,:,:) 
    192       IF( num_sal == 2 )   smv_i(:,:,:) = old_smv_i(:,:,:) 
    193       ! 
    194       IF( nn_timing == 1 )  CALL timing_stop('limitd_th') 
     167     IF( nn_timing == 1 )  CALL timing_stop('limitd_th') 
    195168   END SUBROUTINE lim_itd_th 
    196169   ! 
    … …  
    281254         DO jj = 1, jpj 
    282255            DO ji = 1, jpi 
    283                zindb             = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)+epsi10))     !0 if no ice and 1 if yes 
    284                ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi10) * zindb 
    285                zindb             = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 
    286                zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),epsi10) * zindb 
    287                IF( a_i(ji,jj,jl) > 1e-6 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)  
     256               zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes 
     257               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb 
     258               zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 
     259               zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb 
     260               IF( a_i(ji,jj,jl) > epsi06 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)  
    288261            END DO 
    289262         END DO 
    … …  
    321294      ! Tricky trick see limitd_me.F90 
    322295      ! will be soon removed, CT 
    323       ! hi_max(kubnd) = 999.99 
     296      ! hi_max(kubnd) = 99. 
    324297      zhbnew(:,:,:) = 0._wp 
    325298 
    … …  
    329302            ij = nind_j(ji) 
    330303            ! 
    331             IF ( ( zht_i_o(ii,ij,jl)  .GT.epsi10 ) .AND. &  
    332                ( zht_i_o(ii,ij,jl+1).GT.epsi10 ) ) THEN 
     304            IF ( ( zht_i_o(ii,ij,jl) .GT. epsi10 ) .AND. &  
     305               ( zht_i_o(ii,ij,jl+1) .GT. epsi10 ) ) THEN 
    333306               !interpolate between adjacent category growth rates 
    334307               zslope = ( zdhice(ii,ij,jl+1)     - zdhice(ii,ij,jl) ) / & 
    … …  
    402375               zhbnew(ji,jj,kubnd) = 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) 
    403376            ELSE 
    404                zhbnew(ji,jj,kubnd) = hi_max(kubnd) 
     377               zhbnew(ji,jj,kubnd) = hi_max(kubnd)   
     378               !!? clem bug: since hi_max(jpl)=99, this limit is very high  
     379               !!? but I think it is erased in fitline subroutine  
    405380            ENDIF 
    406381 
    … …  
    447422                     * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 
    448423                  a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0 
    449                   v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem@useless ? 
     424                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ? 
    450425               ENDIF     ! zetamax > 0 
    451426               ! ji, a_i > epsi10 
    … …  
    539514            a_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim  
    540515            ht_i(ii,ij,1) = hiclim 
    541             v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem@useless 
     516            v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem-useless 
    542517         ENDIF 
    543518      END DO !ji 
    … …  
    602577      REAL(wp) ::   zdhr         ! 1 / (hR - hL) 
    603578      REAL(wp) ::   zwk1, zwk2   ! temporary variables 
    604       REAL(wp) ::   zacrith      ! critical minimum concentration in an ice category 
    605       !!------------------------------------------------------------------ 
    606       ! 
    607       zacrith       = 1.0e-6 
     579      !!------------------------------------------------------------------ 
     580      ! 
    608581      ! 
    609582      DO jj = 1, jpj 
    610583         DO ji = 1, jpi 
    611584            ! 
    612             IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > zacrith   & 
     585            IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10   & 
    613586               &                        .AND. hice(ji,jj)        > 0._wp     ) THEN 
    614587 
    … …  
    1009982                  !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) 
    1010983                  !zdvice(ji,jj,jl)  = v_i(ji,jj,jl) 
    1011                   zdaice(ji,jj,jl)  = a_i(ji,jj,jl)/2 
    1012                   zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1))/2 
     984                  !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * 0.5_wp 
     985                  !zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp 
    1013986                  ! end TECLIM change  
     987                  ! clem: how much of a_i you send in cat sup is somewhat arbitrary 
     988                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi10 ) / ht_i(ji,jj,jl)   
     989                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 ) 
    1014990               ENDIF 
    1015991            END DO                 ! ji 
    … …  
    10381014         zshiftflag = 0 
    10391015 
     1016!clem-change 
     1017!         DO jj = 1, jpj 
     1018!            DO ji = 1, jpi 
     1019!               IF( a_i(ji,jj,jl+1) >  epsi10 .AND.   & 
     1020!                  ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
     1021!                  ! 
     1022!                  zshiftflag = 1 
     1023!                  zdonor(ji,jj,jl) = jl + 1 
     1024!                  zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)  
     1025!                  zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 
     1026!               ENDIF 
     1027!            END DO                 ! ji 
     1028!         END DO                 ! jj 
     1029! 
     1030!         IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
     1031!          
     1032!         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
     1033!            CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
     1034!            ! Reset shift parameters 
     1035!            zdonor(:,:,jl) = 0 
     1036!            zdaice(:,:,jl) = 0._wp 
     1037!            zdvice(:,:,jl) = 0._wp 
     1038!         ENDIF 
     1039!clem-change 
     1040 
     1041         ! clem-change begin: why not doing that? 
    10401042         DO jj = 1, jpj 
    10411043            DO ji = 1, jpi 
    10421044               IF( a_i(ji,jj,jl+1) >  epsi10 .AND.   & 
    10431045                  ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
    1044                   ! 
    1045                   zshiftflag = 1 
    1046                   zdonor(ji,jj,jl) = jl + 1 
    1047                   zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)  
    1048                   zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 
     1046                  ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 
     1047                  a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)  
    10491048               ENDIF 
    10501049            END DO                 ! ji 
    10511050         END DO                 ! jj 
    1052  
    1053          IF(lk_mpp)   CALL mpp_max( zshiftflag ) 
    1054           
    1055          IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
    1056             CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
    1057             ! Reset shift parameters 
    1058             zdonor(:,:,jl) = 0 
    1059             zdaice(:,:,jl) = 0._wp 
    1060             zdvice(:,:,jl) = 0._wp 
    1061          ENDIF 
     1051         ! clem-change end 
    10621052 
    10631053      END DO                    ! jl 

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

    r4220 r4332  
    5050   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2) 
    5151 
     52   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    5253   REAL(wp) ::   rzero   = 0._wp   ! constant values 
    5354   REAL(wp) ::   rone    = 1._wp   ! constant values 
    … …  
    486487         CALL lbc_lnk( zs12(:,:), 'F', 1. ) 
    487488 
    488 !#if defined key_bdy 
    489 !         ! clem: change zs1, zs2, zs12 at the boundary for each iteration 
    490 !         CALL bdy_ice_lim_dyn( 2, zs1, zs2, zs12 ) 
    491 !         CALL lbc_lnk( zs1 (:,:), 'T', 1. ) 
    492 !         CALL lbc_lnk( zs2 (:,:), 'T', 1. ) 
    493 !         CALL lbc_lnk( zs12(:,:), 'F', 1. ) 
    494 !#endif          
    495  
    496489         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
    497490         DO jj = k_j1+1, k_jpj-1 
    … …  
    544537            CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
    545538#endif 
     539#if defined key_bdy 
     540         ! clem: change u_ice and v_ice at the boundary for each iteration 
     541         CALL bdy_ice_lim_dyn( 'U' ) 
     542#endif          
    546543 
    547544!CDIR NOVERRCHK 
    … …  
    572569            CALL agrif_rhg_lim2( jter, nevp, 'V' ) 
    573570#endif 
     571#if defined key_bdy 
     572         ! clem: change u_ice and v_ice at the boundary for each iteration 
     573         CALL bdy_ice_lim_dyn( 'V' ) 
     574#endif          
    574575 
    575576         ELSE  
    … …  
    599600            CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    600601#if defined key_agrif && defined key_lim2 
    601             CALL agrif_rhg_lim2( jter, nevp , 'V' ) 
    602 #endif 
     602            CALL agrif_rhg_lim2( jter, nevp, 'V' ) 
     603#endif 
     604#if defined key_bdy 
     605         ! clem: change u_ice and v_ice at the boundary for each iteration 
     606         CALL bdy_ice_lim_dyn( 'V' ) 
     607#endif          
    603608 
    604609!CDIR NOVERRCHK 
    … …  
    632637            CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
    633638#endif 
     639#if defined key_bdy 
     640         ! clem: change u_ice and v_ice at the boundary for each iteration 
     641         CALL bdy_ice_lim_dyn( 'U' ) 
     642#endif          
    634643 
    635644         ENDIF 
    636645          
    637 !#if defined key_bdy 
    638 !         ! clem: change u_ice and v_ice at the boundary for each iteration 
    639 !         CALL bdy_ice_lim_dyn( 1 ) 
    640 !#endif          
    641  
    642646         IF(ln_ctl) THEN 
    643647            !---  Convergence test. 
    … …  
    666670!CDIR NOVERRCHK 
    667671         DO ji = fs_2, fs_jpim1 
    668             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    669             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     672            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
     673            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    670674            zdummy = vt_i(ji,jj) 
    671675            IF ( zdummy .LE. hminrhg ) THEN 
    … …  
    684688#if defined key_bdy 
    685689      ! clem: change u_ice and v_ice at the boundary 
    686       CALL bdy_ice_lim_dyn( 1 ) 
     690      CALL bdy_ice_lim_dyn( 'U' ) 
     691      CALL bdy_ice_lim_dyn( 'V' ) 
    687692#endif          
    688693 
    689694      DO jj = k_j1+1, k_jpj-1  
    690695         DO ji = fs_2, fs_jpim1 
    691             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    692             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     696            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
     697            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    693698            zdummy = vt_i(ji,jj) 
    694699            IF ( zdummy .LE. hminrhg ) THEN 
    … …  
    714719            !- zdd(:,:), zdt(:,:): divergence and tension at centre  
    715720            !- zds(:,:): shear on northeast corner of grid cells 
    716             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    717             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     721            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
     722            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    718723            zdummy = vt_i(ji,jj) 
    719724            IF ( zdummy .LE. hminrhg ) THEN 

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

    r4220 r4332  
    4949   PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim 
    5050 
    51    REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values 
    5251   REAL(wp)  ::   rzero  = 0._wp     
    5352   REAL(wp)  ::   rone   = 1._wp 
    … …  
    145144 
    146145            !   computation the solar flux at ocean surface 
    147             IF (lk_cpl) THEN ! be carfeful: not being tested yet 
     146            IF (lk_cpl) THEN ! be carfeful: not been tested yet 
    148147               ! original line 
    149148               !zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) 
    … …  
    188187            qns(ji,jj) = zfcm2 - fdtcn(ji,jj)                        ! non solar heat flux 
    189188            !                           ! fdtcn : turbulent oceanic heat flux 
    190  
    191             !!gm   this IF prevents the vertorisation of the whole loop 
    192           !  IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 
    193           !     WRITE(numout,*) ' lim_sbc : heat fluxes ' 
    194           !     WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx) 
    195           !     WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx) 
    196           !     WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx) 
    197           !     WRITE(numout,*) 
    198           !     WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx) 
    199           !     WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
    200           !     WRITE(numout,*) ' ifral     : ', ifral 
    201           !     WRITE(numout,*) ' ial       : ', ial   
    202           !     WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx) 
    203           !     WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx) 
    204           !     !WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice 
    205           !     !WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice 
    206           !     WRITE(numout,*) ' ifrdv     : ', ifrdv 
    207           !     WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx) 
    208           !     WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx) 
    209           !     !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice 
    210           !     !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice 
    211           !     WRITE(numout,*) ' ' 
    212           !     WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
    213           !     WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx) 
    214           !     WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx) 
    215           !     WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx) 
    216           !     WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 
    217           !  ENDIF 
    218             !!gm   end 
    219189         END DO 
    220190      END DO 

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

    r4220 r4332  
    5050   PUBLIC   lim_thd_init   ! called by iceini module 
    5151 
    52    REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values 
    53    REAL(wp) ::   epsi16 = 1e-16_wp   ! 
    54    REAL(wp) ::   epsi10 = 1e-10_wp   ! 
    55    REAL(wp) ::   epsi06 = 1e-06_wp   ! 
    56    REAL(wp) ::   epsi04 = 1e-04_wp   ! 
     52   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    5753   REAL(wp) ::   zzero  = 0._wp      ! 
    5854   REAL(wp) ::   zone   = 1._wp      ! 
    … …  
    126122               DO ji = 1, jpi 
    127123                  !Energy of melting q(S,T) [J.m-3] 
    128                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * REAL( nlay_i ) 
     124                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 
    129125                  !0 if no ice and 1 if yes 
    130                   zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) )  )  
     126                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 )  ) 
    131127                  !convert units ! very important that this line is here 
    132128                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb  
    … …  
    138134               DO ji = 1, jpi 
    139135                  !Energy of melting q(S,T) [J.m-3] 
    140                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL( nlay_s ) 
     136                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 
    141137                  !0 if no ice and 1 if yes 
    142                   zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) )  )  
     138                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 )  ) 
    143139                  !convert units ! very important that this line is here 
    144140                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb  
    … …  
    147143         END DO 
    148144      END DO 
    149  
    150       !----------------------------- 
    151       ! 1.3) Set some dummies to 0 
    152       !----------------------------- 
    153       !clem rdvosif(:,:) = 0.e0   ! variation of ice volume at surface 
    154       !clem rdvobif(:,:) = 0.e0   ! variation of ice volume at bottom 
    155       !clem fdvolif(:,:) = 0.e0   ! total variation of ice volume 
    156       !clem rdvonif(:,:) = 0.e0   ! lateral variation of ice volume 
    157       !clem fstric (:,:) = 0.e0   ! part of solar radiation transmitted through the ice 
    158       !clem ffltbif(:,:) = 0.e0   ! linked with fstric 
    159       !clem qfvbq  (:,:) = 0.e0   ! linked with fstric 
    160       !clem rdm_snw(:,:) = 0.e0   ! variation of snow mass per unit area 
    161       !clem rdm_ice(:,:) = 0.e0   ! variation of ice mass per unit area 
    162       !clem hicifp (:,:) = 0.e0   ! daily thermodynamic ice production.  
    163       !clem sfx_bri(:,:) = 0.e0   ! brine flux contribution to salt flux to the ocean 
    164       !clem fhbri  (:,:) = 0.e0   ! brine flux contribution to heat flux to the ocean 
    165       !clem sfx_thd(:,:) = 0.e0   ! equivalent salt flux to the ocean due to ice/growth decay 
    166145 
    167146      !----------------------------------- 
    … …  
    182161!CDIR NOVERRCHK 
    183162         DO ji = 1, jpi 
    184             phicif(ji,jj)  = vt_i(ji,jj) 
    185             pfrld(ji,jj)   = 1.0 - at_i(ji,jj) 
    186             zinda          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) ) 
     163            zinda          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) + epsi10 ) ) ) 
    187164            ! 
    188165            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
    … …  
    195172 
    196173            ! here the drag will depend on ice thickness and type (0.006) 
    197             fdtcn(ji,jj)  = zinda * rau0 * rcp * 0.006  * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) )  
     174            fdtcn(ji,jj)  = zinda * rau0 * rcp * 0.006  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) )  
    198175            ! also category dependent 
    199176            !           !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead  
    200             qdtcn(ji,jj)  = zinda * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 
     177            qdtcn(ji,jj)  = zinda * fdtcn(ji,jj) * ( 1.0 - at_i(ji,jj) ) * rdt_ice 
    201178            !                        
    202179            !           !-- Lead heat budget, qldif (part 1, next one is in limthd_dh)  
    … …  
    214191            zpareff        = 1.0 - zinda * zfntlat 
    215192            != 0 if ice and positive heat budget and 1 if one of those two is false 
    216             zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( at_i(ji,jj) * rdt_ice , epsi16 ) 
     193            zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / ( rdt_ice * MAX( at_i(ji,jj), epsi10 ) ) 
    217194            ! 
    218195            ! Heat budget of the lead, energy transferred from ice to ocean 
    … …  
    221198            ! 
    222199            ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 
    223             qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) 
     200            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 
    224201            ! 
    225202            ! oceanic heat flux (limthd_dh) 
    226             fbif   (ji,jj) = zinda * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 
     203            fbif   (ji,jj) = zinda * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) + fdtcn(ji,jj) ) 
    227204            ! 
    228205         END DO 
    … …  
    240217         ENDIF 
    241218 
    242          zareamin = 1.e-10 
     219         zareamin = epsi10 
    243220         nbpb = 0 
    244221         DO jj = 1, jpj 
    … …  
    248225                  npb(nbpb) = (jj - 1) * jpi + ji 
    249226               ENDIF 
    250                ! debug point to follow 
    251                IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN 
    252                   jiindex_1d = nbpb 
    253                ENDIF 
    254             END DO 
    255          END DO 
     227            END DO 
     228         END DO 
     229 
     230         ! debug point to follow 
     231         jiindex_1d = 0 
     232         IF( ln_nicep ) THEN 
     233            DO ji = mi0(jiindx), mi1(jiindx) 
     234               DO jj = mj0(jjindx), mj1(jjindx) 
     235                  jiindex_1d = (jj - 1) * jpi + ji 
     236               END DO 
     237            END DO 
     238         ENDIF 
    256239 
    257240         !------------------------------------------------------------------------------! 
    … …  
    315298            !-------------------------------- 
    316299 
    317             IF( con_i )   CALL lim_thd_enmelt( 1, nbpb )   ! computes sea ice energy of melting 
    318             IF( con_i )   CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl ) 
     300            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_enmelt( 1, nbpb )   ! computes sea ice energy of melting 
     301            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_glohec( qt_i_in, qt_s_in, q_i_layer_in, 1, nbpb, jl ) 
    319302 
    320303            !                                 !---------------------------------! 
    … …  
    324307            CALL lim_thd_enmelt( 1, nbpb )    ! computes sea ice energy of melting compulsory for limthd_dh 
    325308 
    326             IF( con_i )   CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )  
    327             IF( con_i )   CALL lim_thd_con_dif( 1 , nbpb , jl ) 
     309            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_glohec ( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )  
     310            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_con_dif( 1 , nbpb , jl ) 
    328311 
    329312            !                                 !---------------------------------! 
    … …  
    340323 
    341324            !           CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting 
    342             IF( con_i )   CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )  
    343             IF( con_i )   CALL lim_thd_con_dh ( 1 , nbpb , jl ) 
     325            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_glohec( qt_i_fin, qt_s_fin, q_i_layer_fin, 1, nbpb, jl )  
     326            IF( con_i .AND. jiindex_1d > 0 )   CALL lim_thd_con_dh ( 1 , nbpb , jl ) 
    344327 
    345328            !-------------------------------- 
    … …  
    431414!clem@mv-to-itd    dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    432415 
    433       IF( con_i )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 
     416      IF( con_i .AND. jiindex_1d > 0 )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 
    434417 
    435418      IF(ln_ctl) THEN            ! Control print 
    … …  
    522505      END DO 
    523506      ! 
    524       IF(lwp) WRITE(numout,*) ' lim_thd_glohec ' 
    525       IF(lwp) WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice 
    526       IF(lwp) WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice 
    527       IF(lwp) WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice 
     507      WRITE(numout,*) ' lim_thd_glohec ' 
     508      WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice 
     509      WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice 
     510      WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice 
    528511      ! 
    529512   END SUBROUTINE lim_thd_glohec 
    … …  
    611594      WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit 
    612595 
    613       IF (jiindex_1D.GT.0) WRITE(numout,*) ' fc_su      : ', fc_su(jiindex_1d) 
    614       IF (jiindex_1D.GT.0) WRITE(numout,*) ' fatm       : ', fatm(jiindex_1d,jl) 
    615       IF (jiindex_1D.GT.0) WRITE(numout,*) ' t_su       : ', t_su_b(jiindex_1d) 
     596      WRITE(numout,*) ' fc_su      : ', fc_su(jiindex_1d) 
     597      WRITE(numout,*) ' fatm       : ', fatm(jiindex_1d,jl) 
     598      WRITE(numout,*) ' t_su       : ', t_su_b(jiindex_1d) 
    616599 
    617600      !--------------------------------------- 

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

    r4099 r4332  
    3232   PUBLIC   lim_thd_dh   ! called by lim_thd 
    3333 
    34    REAL(wp) ::   epsi20 = 1e-20   ! constant values 
    35    REAL(wp) ::   epsi13 = 1e-13   ! 
    36    REAL(wp) ::   epsi16 = 1e-16   ! 
    37    REAL(wp) ::   zzero  = 0.e0    ! 
    38    REAL(wp) ::   zone   = 1.e0    ! 
     34   REAL(wp) ::   epsi20 = 1.e-20   ! constant values 
     35   REAL(wp) ::   epsi10 = 1.e-10   ! 
     36   REAL(wp) ::   epsi13 = 1.e-13   ! 
     37   REAL(wp) ::   zzero  = 0._wp    ! 
     38   REAL(wp) ::   zone   = 1._wp    ! 
    3939 
    4040   !!---------------------------------------------------------------------- 
    … …  
    125125      INTEGER  ::   num_iter_max, numce_dh 
    126126      REAL(wp) ::   meance_dh 
     127      REAL(wp) ::   zinda  
    127128      REAL(wp), POINTER, DIMENSION(:) ::   zinnermelt 
    128129      REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i 
    … …  
    288289      END DO 
    289290 
    290       !                     !------------------- 
    291       IF( con_i ) THEN      ! Conservation test 
    292          !                  !------------------- 
     291      !                                          !------------------- 
     292      IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test 
     293         !                                       !------------------- 
    293294         numce_dh  = 0 
    294295         meance_dh = 0._wp 
    … …  
    494495      END DO ! jk 
    495496 
    496       !                     !------------------- 
    497       IF( con_i ) THEN      ! Conservation test 
    498       !                     !------------------- 
     497      !                                          !------------------- 
     498      IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test 
     499      !                                          !------------------- 
    499500         DO ji = kideb, kiut 
    500501            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0  ) THEN 
    … …  
    544545         sfx_thd_1d(ji)  = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 
    545546         !                     ! excessive energy is sent to lateral ablation 
    546          fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres * r1_rdtice 
     547         zinda = MAX( 0._wp, SIGN( 1._wp , 1.0 - at_i_b(ji) - epsi10 ) ) 
     548         fsup(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi10 ) * zdvres * r1_rdtice 
    547549      END DO 
    548550 
    … …  
    577579 
    578580         ! energy of melting of remaining snow 
    579          zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 ) 
     581         zinda = MAX( 0._wp, SIGN( 1._wp , zhni - epsi10 ) ) 
     582         zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi10 ) * zinda 
    580583         zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 ) 
    581584         zhnfi     =  zhni + zdhnm 
    … …  
    613616         ENDIF 
    614617         ! 
    615          ! diagnostic ( bottom ice growth ) 
     618         ! diagnostic  
    616619         ii = MOD( npb(ji) - 1, jpi ) + 1 
    617620         ij = ( npb(ji) - 1 ) / jpi + 1 
    … …  
    689692         ! clem: new salinity difference stored (to be used in limthd_ent.F90) 
    690693         IF (  num_sal == 2  ) THEN 
    691             i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) ) 
     694            i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) ) 
    692695            ! salinity dif due to snow-ice formation 
    693             dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch      
     696            dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch      
    694697            ! salinity dif due to bottom growth  
    695698            IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN 
    696                dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 
     699               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch 
    697700            ENDIF 
    698701         ENDIF 

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

    r4220 r4332  
    3131   PUBLIC   lim_thd_dif   ! called by lim_thd 
    3232 
    33    REAL(wp) ::   epsi20 = 1e-20     ! constant values 
    34    REAL(wp) ::   epsi13 = 1e-13     ! constant values 
    35  
     33   REAL(wp) ::   epsi10      =  1.e-10_wp    ! 
    3634   !!---------------------------------------------------------------------- 
    3735   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    … …  
    107105      INTEGER, DIMENSION(kiut) ::   numeqmax   ! reference number of bottom equation 
    108106      INTEGER, DIMENSION(kiut) ::   isnow      ! switch for presence (1) or absence (0) of snow 
    109       REAL(wp) ::   zeps      =  1.e-10_wp    ! 
    110107      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    111108      REAL(wp) ::   zg1       =  2._wp        ! 
    … …  
    114111      REAL(wp) ::   zraext_s  =  1.e+8_wp     ! extinction coefficient of radiation in the snow 
    115112      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity 
    116       REAL(wp) ::   zht_smin  =  1.e-4_wp     ! minimum snow depth 
    117113      REAL(wp) ::   ztmelt_i    ! ice melting temperature 
    118114      REAL(wp) ::   zerritmax   ! current maximal error on temperature  
    … …  
    312308         IF( thcon_i_swi == 0 ) THEN      ! Untersteiner (1964) formula 
    313309            DO ji = kideb , kiut 
    314                ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-zeps,t_i_b(ji,1)-rtt) 
     310               ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt) 
    315311               ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin) 
    316312            END DO 
    … …  
    318314               DO ji = kideb , kiut 
    319315                  ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) /  & 
    320                      MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) 
     316                     MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) 
    321317                  ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) 
    322318               END DO 
    … …  
    326322         IF( thcon_i_swi == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 
    327323            DO ji = kideb , kiut 
    328                ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -zeps, t_i_b(ji,1)-rtt )   & 
     324               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt )   & 
    329325                  &                   - 0.011_wp * ( t_i_b(ji,1) - rtt )   
    330326               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
    … …  
    333329               DO ji = kideb , kiut 
    334330                  ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   & 
    335                      &                                  / MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   & 
     331                     &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   & 
    336332                     &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )   
    337333                  ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 
    … …  
    339335            END DO 
    340336            DO ji = kideb , kiut 
    341                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-zeps,t_bo_b(ji)-rtt)   & 
     337               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   & 
    342338                  &                        - 0.011_wp * ( t_bo_b(ji) - rtt )   
    343339               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 
    … …  
    352348 
    353349            !-- Snow kappa factors 
    354             zkappa_s(ji,0)         = rcdsn / MAX(zeps,zh_s(ji)) 
    355             zkappa_s(ji,nlay_s)    = rcdsn / MAX(zeps,zh_s(ji)) 
     350            zkappa_s(ji,0)         = rcdsn / MAX(epsi10,zh_s(ji)) 
     351            zkappa_s(ji,nlay_s)    = rcdsn / MAX(epsi10,zh_s(ji)) 
    356352         END DO 
    357353 
    … …  
    359355            DO ji = kideb , kiut 
    360356               zkappa_s(ji,layer)  = 2.0 * rcdsn / & 
    361                   MAX(zeps,2.0*zh_s(ji)) 
     357                  MAX(epsi10,2.0*zh_s(ji)) 
    362358            END DO 
    363359         END DO 
    … …  
    367363               !-- Ice kappa factors 
    368364               zkappa_i(ji,layer)  = 2.0*ztcond_i(ji,layer)/ & 
    369                   MAX(zeps,2.0*zh_i(ji))  
    370             END DO 
    371          END DO 
    372  
    373          DO ji = kideb , kiut 
    374             zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(zeps,zh_i(ji)) 
    375             zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji)) 
     365                  MAX(epsi10,2.0*zh_i(ji))  
     366            END DO 
     367         END DO 
     368 
     369         DO ji = kideb , kiut 
     370            zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji)) 
     371            zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji)) 
    376372            !-- Interface 
    377             zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & 
     373            zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, & 
    378374               (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 
    379375            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 
    … …  
    389385               ztitemp(ji,layer)   = t_i_b(ji,layer) 
    390386               zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ & 
    391                   MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),zeps) 
     387                  MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10) 
    392388               zeta_i(ji,layer)    = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), & 
    393                   zeps) 
     389                  epsi10) 
    394390            END DO 
    395391         END DO 
    … …  
    398394            DO ji = kideb , kiut 
    399395               ztstemp(ji,layer) = t_s_b(ji,layer) 
    400                zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),zeps) 
     396               zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 
    401397            END DO 
    402398         END DO 
    … …  
    707703      END DO  ! End of the do while iterative procedure 
    708704 
    709       IF( ln_nicep ) THEN 
     705      IF( ln_nicep .AND. lwp ) THEN 
    710706         WRITE(numout,*) ' zerritmax : ', zerritmax 
    711707         WRITE(numout,*) ' nconv     : ', nconv 
    … …  
    732728      ! Heat conservation       ! 
    733729      !-------------------------! 
    734       IF( con_i ) THEN 
     730      IF( con_i .AND. jiindex_1d > 0 ) THEN 
    735731         DO ji = kideb, kiut 
    736732            ! Upper snow value 

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

    r4072 r4332  
    3636   PUBLIC   lim_thd_ent         ! called by lim_thd 
    3737 
    38    REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values 
    39    REAL(wp) ::   epsi13 = 1e-13_wp   ! 
    40    REAL(wp) ::   epsi10 = 1e-10_wp   ! 
    41    REAL(wp) ::   epsi06 = 1e-06_wp   ! 
     38   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values 
     39   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    4240   REAL(wp) ::   zzero  = 0._wp      ! 
    4341   REAL(wp) ::   zone   = 1._wp      ! 
    … …  
    122120      REAL(wp), POINTER, DIMENSION(:,:) ::   zhl0      ! old and new layer thicknesses  
    123121      REAL(wp), POINTER, DIMENSION(:,:) ::   zrl01 
     122 
     123      REAL(wp) ::   zinda  
    124124      !!------------------------------------------------------------------- 
    125125 
    … …  
    391391         DO layer1 = ntop1, nbot1 
    392392            DO ji = kideb, kiut 
    393                zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1))   & 
     393               zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) 
     394               zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1))   & 
    394395                  &                 - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10))  
    395396               q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0)   & 
    … …  
    407408      END DO 
    408409 
    409       IF ( con_i ) THEN 
     410      IF ( con_i .AND. jiindex_1d > 0 ) THEN 
    410411         DO ji = kideb, kiut 
    411412            IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN 
    … …  
    429430      DO jk = 1, nlay_s 
    430431         DO ji = kideb, kiut 
    431             q_s_b(ji,jk) = q_s_b(ji,jk) / MAX( zh_s(ji) , epsi20 ) 
     432            zinda = MAX( 0._wp, SIGN( 1._wp , zh_s(ji) - epsi10 ) )         
     433            q_s_b(ji,jk) = zinda * q_s_b(ji,jk) / MAX( zh_s(ji) , epsi10 ) 
    432434         END DO !ji 
    433435      END DO !jk   
    … …  
    555557         DO ji = kideb, kiut 
    556558            ! Heat conservation 
    557             zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06) ) & 
     559            zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi10) ) & 
    558560               &                                   * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) 
    559561         END DO 
    … …  
    601603         DO layer1 = ntop1, nbot1 
    602604            DO ji = kideb, kiut 
    603                zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) & 
     605               zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) 
     606               zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) & 
    604607                  - MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10)) 
    605608               q_i_b(ji,layer1) = q_i_b(ji,layer1) &  
    606609                  + zrl01(layer1,layer0)*qm0(ji,layer0) & 
    607                   * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06)) & 
     610                  * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi10)) & 
    608611                  * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 
    609612            END DO 
    … …  
    621624      END DO 
    622625      ! 
    623       IF ( con_i ) THEN 
     626      IF ( con_i .AND. jiindex_1d > 0 ) THEN 
    624627         DO ji = kideb, kiut 
    625628            IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN 
    … …  
    646649      DO jk = 1, nlay_i 
    647650         DO ji = kideb, kiut 
    648             q_i_b(ji,jk) = q_i_b(ji,jk) / MAX( zh_i(ji) , epsi20 ) 
     651            zinda = MAX( 0._wp, SIGN( 1._wp , zh_i(ji) - epsi10 ) ) 
     652            q_i_b(ji,jk) = zinda * q_i_b(ji,jk) / MAX( zh_i(ji) , epsi10 ) 
    649653         END DO !ji 
    650654      END DO !jk   

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

    r4045 r4332  
    3636   PUBLIC lim_thd_lac     ! called by lim_thd 
    3737 
    38    REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values 
    39    REAL(wp) ::   epsi13 = 1e-13_wp   ! 
    40    REAL(wp) ::   epsi11 = 1e-11_wp   ! 
    41    REAL(wp) ::   epsi10 = 1e-10_wp   ! 
    42    REAL(wp) ::   epsi06 = 1e-06_wp   ! 
    43    REAL(wp) ::   epsi03 = 1e-03_wp   ! 
     38   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    4439   REAL(wp) ::   zzero  = 0._wp      ! 
    4540   REAL(wp) ::   zone   = 1._wp      ! 
    … …  
    160155                  !Energy of melting q(S,T) [J.m-3] 
    161156                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * REAL( nlay_i ) 
    162                   zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) )  )   !0 if no ice and 1 if yes 
     157                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes 
    163158                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 
    164159               END DO 
    … …  
    286281               nbpac = nbpac + 1 
    287282               npac( nbpac ) = (jj - 1) * jpi + ji 
    288                IF( ji == jiindx .AND. jj == jjindx )   jiindex_1d = nbpac 
    289283            ENDIF 
    290284         END DO 
    291285      END DO 
    292286 
    293       IF( ln_nicep )   WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 
     287      ! debug point to follow 
     288      jiindex_1d = 0 
     289      IF( ln_nicep ) THEN 
     290         DO ji = mi0(jiindx), mi1(jiindx) 
     291            DO jj = mj0(jjindx), mj1(jjindx) 
     292               IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN 
     293                  jiindex_1d = (jj - 1) * jpi + ji 
     294               ENDIF 
     295            END DO 
     296         END DO 
     297      ENDIF 
     298    
     299      IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 
    294300 
    295301      !------------------------------ 
    … …  
    499505         END DO 
    500506 
    501          IF( ln_nicep ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl) 
     507         IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 
    502508         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    503509            DO ji = 1, nbpac 
    504510               zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) ) 
    505                zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi06 ) )  ! clem 
    506                zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 ) 
    507             END DO 
    508          END DO 
    509          IF( ln_nicep )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl) 
     511               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) )  ! clem 
     512               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 ) 
     513            END DO 
     514         END DO 
     515         IF( ln_nicep .AND. jiindex_1d > 0 )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 
    510516 
    511517         !--------------------------------- 
    … …  
    649655         !     CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)  
    650656         IF( ln_nicep ) THEN 
    651             WRITE(numout,*) ' vt_i_init : ', vt_i_init(jiindx,jjindx) 
    652             WRITE(numout,*) ' vt_i_final: ', vt_i_final(jiindx,jjindx) 
    653             WRITE(numout,*) ' et_i_init : ', et_i_init(jiindx,jjindx) 
    654             WRITE(numout,*) ' et_i_final: ', et_i_final(jiindx,jjindx) 
     657            DO ji = mi0(jiindx), mi1(jiindx) 
     658               DO jj = mj0(jjindx), mj1(jjindx) 
     659                  WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj) 
     660                  WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj) 
     661                  WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj) 
     662                  WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj) 
     663               END DO 
     664            END DO 
    655665         ENDIF 
    656666         ! 

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

    r4220 r4332  
    3636   PUBLIC   lim_trp    ! called by ice_step 
    3737 
    38    REAL(wp)  ::   epsi06 = 1.e-06_wp   ! constant values 
    39    REAL(wp)  ::   epsi03 = 1.e-03_wp   
    4038   REAL(wp)  ::   epsi10 = 1.e-10_wp   
    41    REAL(wp)  ::   epsi16 = 1.e-16_wp 
    42    REAL(wp)  ::   epsi20 = 1.e-20_wp 
    4339   REAL(wp)  ::   rzero  = 0._wp    
    4440   REAL(wp)  ::   rone   = 1._wp 
    45  
    46    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:,:) ::   zs0e 
    4741 
    4842   !! * Substitution 
    … …  
    448442 
    449443                  ! Switches and dummy variables 
    450                   zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi16 ) 
    451                   zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi16 ) 
     444                  zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi10 ) 
     445                  zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi10 ) 
    452446                  zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 
    453447                  zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 
    … …  
    460454                  ENDIF 
    461455 
    462                   zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi16 ) ), 0._wp  ) * a_i(ji,jj,jl) 
     456                  zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) ), 0._wp  ) * a_i(ji,jj,jl) 
    463457                  oa_i (ji,jj,jl)  = zindic * zage  
    464458 
    … …  
    500494               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 
    501495               ! 
    502                zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi16 ) ) 
    503                icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda  ! ice thickness 
     496               zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi10 ) ) 
     497               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda  ! ice thickness 
    504498            END DO 
    505499         END DO 

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

    r4099 r4332  
    4848   PUBLIC   lim_update1   ! routine called by ice_step 
    4949 
    50       REAL(wp)  ::   epsi06 = 1.e-06_wp   ! module constants 
    51       REAL(wp)  ::   epsi04 = 1.e-04_wp   !    -       - 
    52       REAL(wp)  ::   epsi03 = 1.e-03_wp   !    -       - 
    5350      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       - 
    54       REAL(wp)  ::   epsi16 = 1.e-16_wp   !    -       - 
    55       REAL(wp)  ::   epsi20 = 1.e-20_wp   !    -       - 
    5651      REAL(wp)  ::   rzero  = 0._wp       !    -       - 
    5752      REAL(wp)  ::   rone   = 1._wp       !    -       - 
    … …  
    108103      !------------------------------------------------------------------------------ 
    109104 
    110       !--------------------- 
    111       ! Ice dynamics   
    112       !--------------------- 
    113       u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:) 
    114       v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:) 
    115   
    116       !----------------------------- 
    117       ! Update ice and snow volumes   
    118       !----------------------------- 
    119       DO jl = 1, jpl 
    120          v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_trp(:,:,jl)  
    121          v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_trp(:,:,jl) 
    122       END DO 
    123  
    124  
    125       !--------------------------------------------- 
    126       ! Ice concentration and ice heat content 
    127       !--------------------------------------------- 
    128       a_i (:,:,:)  = a_i (:,:,:) + d_a_i_trp(:,:,:) 
    129       e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:)   
    130  
    131       !------------------------------ 
    132       ! Snow temperature and ice age 
    133       !------------------------------ 
    134       e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_trp (:,:,:,:) 
    135       oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_trp(:,:,:)   
    136  
    137       !-------------- 
    138       ! Ice salinity     
    139       !-------------- 
    140       IF(  num_sal == 2  ) THEN      ! general case 
    141          smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_trp(:,:,:) 
    142       ENDIF 
     105      !----------------- 
     106      !  Trend terms 
     107      !----------------- 
     108      d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:) 
     109      d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:) 
     110      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:) 
     111      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)   
     112      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)    
     113      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)   
     114      d_e_i_trp  (:,:,:,:) = e_i  (:,:,:,:) - old_e_i  (:,:,:,:) 
     115      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
     116      d_smv_i_trp(:,:,:)   = 0._wp 
     117      IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
    143118 
    144119      ! mass and salt flux init (clem) 
    … …  
    160135      CALL lim_var_glo2eqv 
    161136 
    162       !--------------------------------- 
    163       ! Classify the pathological cases 
    164       !--------------------------------- 
    165       ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case) 
    166       ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation) 
    167       ! (5) v_i (old) = 0; d_v_i_trp > 0 (advection of ice in a free-cell) 
    168       ! 
    169       DO jl = 1, jpl 
    170          DO jj = 1, jpj 
    171             DO ji = 1, jpi 
    172                patho_case(ji,jj,jl) = 1 
    173                IF( v_i(ji,jj,jl) .LT. 0.0 ) THEN 
    174                   patho_case(ji,jj,jl) = 3 
    175                ENDIF 
    176                IF( ( old_v_i(ji,jj,jl) .LE. epsi10 ) .AND. ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN 
    177                   patho_case(ji,jj,jl) = 5 ! advection of ice in an ice-free cell 
    178                ENDIF 
    179             END DO 
    180          END DO 
    181       END DO 
    182  
    183137      !-------------------------------------- 
    184138      ! 2. Review of all pathological cases 
    185139      !-------------------------------------- 
    186140 
     141! clem: useless now 
    187142      !------------------------------------------- 
    188143      ! 2.1) Advection of ice in an ice-free cell 
    189144      !------------------------------------------- 
    190145      ! should be removed since it is treated after dynamics now 
    191  
    192 !      !IF( ln_nicep ) THEN   
    193 !         WRITE(numout,*) ' limupdate1 - before h correction ' 
    194 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    195 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    196 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    197 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    198 !      !ENDIF 
     146!      zhimax = 5._wp 
     147!      ! first category 
     148!      DO jj = 1, jpj 
     149!         DO ji = 1, jpi 
     150!            !--- the thickness of such an ice is often out of bounds 
     151!            !--- thus we recompute a new area while conserving ice volume 
     152!            zat_i_old = SUM( old_a_i(ji,jj,:) ) 
     153!            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1) ) - epsi10 ) )  
     154!            IF( ( ABS( d_v_i_trp(ji,jj,1) ) / MAX( ABS( d_a_i_trp(ji,jj,1) ), epsi10 ) * zindb .GT. zhimax ) & 
     155!              &   .AND.( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 
     156!              &   .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
     157!               ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 
     158!               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
     159!            ENDIF 
     160!         END DO 
     161!      END DO 
    199162! 
    200       zhimax = .3_wp 
    201       ! first category 
    202       DO jj = 1, jpj 
    203          DO ji = 1, jpi 
    204             !--- the thickness of such an ice is often out of bounds 
    205             !--- thus we recompute a new area while conserving ice volume 
    206             zat_i_old = SUM( old_a_i(ji,jj,:) ) 
    207             zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1)) - epsi10 ) )  
    208             IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb .GT. zhimax) & 
    209                  .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) & 
    210                  .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
    211                ht_i(ji,jj,1) = hi_max(1) / 2.0 
    212                a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
    213             ENDIF 
    214          END DO 
    215       END DO 
    216  
    217  
    218 !      !IF( ln_nicep ) THEN   
    219 !         at_i(:,:) = 0._wp 
    220 !         DO jl = 1, jpl 
    221 !            at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    222 !         END DO 
    223 !         WRITE(numout,*) ' limupdate1 - after h correction 1 ' 
    224 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    225 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    226 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    227 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    228 !      !ENDIF 
    229  
    230       zhimax = 1._wp 
    231       ! other categories 
    232       DO jl = 2, jpl 
    233          jm = ice_types(jl) 
    234          DO jj = 1, jpj 
    235             DO ji = 1, jpi 
    236                zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) )  
    237                ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
    238                ! it makes problems when the advected volume and concentration do not seem to be  
    239                ! related with each other 
    240                ! the new thickness is sometimes very big! 
    241                ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
    242                ! which of course is plausible 
    243                ! but fuck! it fucks everything up :) 
    244                IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb .GT. zhimax) & 
    245                     .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN 
    246                   ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0 
    247                   a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
    248                ENDIF 
    249             END DO ! ji 
    250          END DO !jj 
    251       END DO !jl 
     163!      zhimax = 20._wp 
     164!      ! other categories 
     165!      DO jl = 2, jpl 
     166!         jm = ice_types(jl) 
     167!         DO jj = 1, jpj 
     168!            DO ji = 1, jpi 
     169!               zindb =  MAX( rzero, SIGN( rone, ABS( d_a_i_trp(ji,jj,jl) ) - epsi10 ) )  
     170!               ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
     171!               ! it makes problems when the advected volume and concentration do not seem to be  
     172!               ! related with each other 
     173!               ! the new thickness is sometimes very big! 
     174!               ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
     175!               ! which of course is plausible 
     176!               ! but fuck! it fucks everything up :) 
     177!               IF ( ( ABS( d_v_i_trp(ji,jj,jl) ) / MAX( ABS( d_a_i_trp(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 
     178!                  &  .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 
     179!                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp 
     180!                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
     181!               ENDIF 
     182!            END DO ! ji 
     183!         END DO !jj 
     184!      END DO !jl 
    252185      
    253186      at_i(:,:) = 0._wp 
    … …  
    255188         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    256189      END DO 
    257  
    258 !      !IF( ln_nicep ) THEN   
    259 !         WRITE(numout,*) ' limupdate1 - after h correction 2 ' 
    260 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    261 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    262 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    263 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    264 !      !ENDIF 
    265190 
    266191      !---------------------------------------------------- 
    … …  
    285210 
    286211               !switches 
    287                zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
     212               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
    288213               !switch = 1 if a_i > 1e-06 and 0 if not 
    289                zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not 
    290                zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not 
     214               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 
     215               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 
    291216               ! bug fix 25 avril 2007 
    292217               zindb         = zindb*zindic 
    … …  
    332257               !-------------------------------------------- 
    333258               ! if greater than 0, ice concentration cannot be smaller than 1e-10 
    334                !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 ) 
    335259               a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 
    336260 
    … …  
    352276            DO jj = 1, jpj  
    353277               DO ji = 1, jpi 
    354                   zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) )  
     278                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )  
    355279                  e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 
    356280               END DO ! ji 
    … …  
    370294         DO ji = 1, jpi 
    371295            z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )         
    372             zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) )  
     296            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) )  
    373297            DO jl  = 1, jpl 
    374                z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb 
     298               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 
    375299               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 
    376300               ! 
    377                zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
    378                ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda 
     301               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
     302               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda 
    379303               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 
    380304            END DO 
    … …  
    412336                     !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) ) 
    413337                     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) 
     338                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 
     339                     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) 
    416340                  END DO ! ji 
    417341               END DO ! jj 
    … …  
    512436            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update1  : old_oa_i    : ') 
    513437            CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ') 
    514             CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update1  : Path. case  : ') 
    515438 
    516439            DO jk = 1, nlay_i 

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

    r4099 r4332  
    4545   PUBLIC   lim_update2   ! routine called by ice_step 
    4646 
    47       REAL(wp)  ::   epsi06 = 1.e-06_wp   ! module constants 
    48       REAL(wp)  ::   epsi04 = 1.e-04_wp   !    -       - 
    4947      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       - 
    50       REAL(wp)  ::   epsi16 = 1.e-16_wp   !    -       - 
    51       REAL(wp)  ::   epsi20 = 1.e-20_wp   !    -       - 
    5248      REAL(wp)  ::   rzero  = 0._wp       !    -       - 
    5349      REAL(wp)  ::   rone   = 1._wp       !    -       - 
    … …  
    10298      CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold )   ! clem 
    10399 
    104       !------------------------------------------------------------------------------ 
    105       ! 1. Update of Global variables                                               | 
    106       !------------------------------------------------------------------------------ 
    107  
    108       !----------------------------- 
    109       ! Update ice and snow volumes   
    110       !----------------------------- 
    111       DO jl = 1, jpl 
    112          v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_thd(:,:,jl)  
    113          v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_thd(:,:,jl) 
    114       END DO 
    115   
    116       !--------------------------------------------- 
    117       ! Ice concentration and ice heat content 
    118       !--------------------------------------------- 
    119       a_i (:,:,:) = a_i (:,:,:) + d_a_i_thd(:,:,:) 
    120       e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:)   
    121   
    122       !------------------------------ 
    123       ! Snow temperature and ice age 
    124       !------------------------------ 
    125       e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_thd (:,:,:,:) 
    126       oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_thd(:,:,:) 
    127  
    128       !-------------- 
    129       ! Ice salinity     
    130       !-------------- 
    131       IF(  num_sal == 2  ) THEN      ! general case 
    132          smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:) 
    133       ENDIF 
     100      !---------------------------------------------------------------------------------------- 
     101      ! 1. Computation of trend terms       
     102      !---------------------------------------------------------------------------------------- 
     103      !- Trend terms 
     104      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:)  
     105      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:) 
     106      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)   
     107      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)  
     108      d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 
     109      !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - old_oa_i (:,:,:) 
     110      d_smv_i_thd(:,:,:) = 0._wp 
     111      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     112      ! diag only (clem) 
     113      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    134114 
    135115      ! mass and salt flux init (clem) 
    … …  
    151131      CALL lim_var_glo2eqv 
    152132 
    153       !--------------------------------- 
    154       ! Classify the pathological cases 
    155       !--------------------------------- 
    156       ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case) 
    157       ! (2) v_i (new) > 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation) 
    158       ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation) 
    159       ! (4) v_i (new) < 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation  
    160       ! with negative advection, very pathological ) 
    161       ! 
    162       DO jl = 1, jpl 
    163          DO jj = 1, jpj 
    164             DO ji = 1, jpi 
    165                patho_case(ji,jj,jl) = 1 
    166                IF( v_i(ji,jj,jl) .GE. 0.0 ) THEN 
    167                   IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN  
    168                      patho_case(ji,jj,jl) = 2 
    169                   ENDIF 
    170                ELSE 
    171                   patho_case(ji,jj,jl) = 3 
    172                   IF( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN  
    173                      patho_case(ji,jj,jl) = 4 
    174                   ENDIF 
    175                ENDIF 
    176              END DO 
    177          END DO 
    178       END DO 
    179  
    180133      !-------------------------------------- 
    181134      ! 2. Review of all pathological cases 
    182135      !-------------------------------------- 
     136 
     137! clem: useless now 
    183138      !------------------------------------------- 
    184139      ! 2.1) Advection of ice in an ice-free cell 
    185140      !------------------------------------------- 
    186141      ! should be removed since it is treated after dynamics now 
    187  
    188 !      !IF( ln_nicep ) THEN   
    189 !         WRITE(numout,*) ' limupdate2 - before h correction ' 
    190 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    191 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    192 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    193 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    194 !      !ENDIF 
    195  
    196       zhimax = .3_wp 
    197       ! first category 
    198       DO jj = 1, jpj 
    199          DO ji = 1, jpi 
    200             !--- the thickness of such an ice is often out of bounds 
    201             !--- thus we recompute a new area while conserving ice volume 
    202             zat_i_old = SUM( old_a_i(ji,jj,:) ) 
    203             zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1)) - epsi10 ) )  
    204             IF ( ( ABS(d_v_i_thd(ji,jj,1))/MAX(ABS(d_a_i_thd(ji,jj,1)),epsi10)*zindb .GT. zhimax) & 
    205                  .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) & 
    206                  .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
    207                ht_i(ji,jj,1) = hi_max(1) / 2.0 
    208                a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
    209             ENDIF 
    210          END DO 
    211       END DO 
    212  
    213 !      !IF( ln_nicep ) THEN   
    214 !         at_i(:,:) = 0._wp 
    215 !         DO jl = 1, jpl 
    216 !            at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     142!      zhimax = 5._wp 
     143!      ! first category 
     144!      DO jj = 1, jpj 
     145!         DO ji = 1, jpi 
     146!            !--- the thickness of such an ice is often out of bounds 
     147!            !--- thus we recompute a new area while conserving ice volume 
     148!            zat_i_old = SUM( old_a_i(ji,jj,:) ) 
     149!            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1) ) - epsi10 ) )  
     150!            IF ( ( ABS( d_v_i_thd(ji,jj,1) ) / MAX( ABS( d_a_i_thd(ji,jj,1) ),epsi10 ) * zindb .GT. zhimax ) & 
     151!               &  .AND. ( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 
     152!               &  .AND. ( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
     153!               ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 
     154!               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
     155!            ENDIF 
    217156!         END DO 
    218 !         WRITE(numout,*) ' limupdate2 - after h correction 1 ' 
    219 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    220 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    221 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    222 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    223 !      !ENDIF 
    224 !  
    225       zhimax = 1._wp 
    226       ! other categories 
    227       DO jl = 2, jpl 
    228          jm = ice_types(jl) 
    229          DO jj = 1, jpj 
    230             DO ji = 1, jpi 
    231                zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_thd(ji,jj,jl)) - epsi10 ) )  
    232               ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
    233                ! it makes problems when the advected volume and concentration do not seem to be  
    234                ! related with each other 
    235                ! the new thickness is sometimes very big! 
    236                ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
    237                ! which of course is plausible 
    238                ! but fuck! it fucks everything up :) 
    239                IF ( (ABS(d_v_i_thd(ji,jj,jl))/MAX(ABS(d_a_i_thd(ji,jj,jl)),epsi10)*zindb .GT. zhimax) & 
    240                     .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN 
    241                   ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0 
    242                   a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
    243                ENDIF 
    244             END DO ! ji 
    245          END DO !jj 
    246       END DO !jl 
     157!      END DO 
     158 
     159!      zhimax = 20._wp 
     160!      ! other categories 
     161!      DO jl = 2, jpl 
     162!         jm = ice_types(jl) 
     163!         DO jj = 1, jpj 
     164!            DO ji = 1, jpi 
     165!               zindb =  MAX( rzero, SIGN( rone, ABS( d_a_i_thd(ji,jj,jl)) - epsi10 ) )  
     166!              ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
     167!               ! it makes problems when the advected volume and concentration do not seem to be  
     168!               ! related with each other 
     169!               ! the new thickness is sometimes very big! 
     170!               ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
     171!               ! which of course is plausible 
     172!               ! but fuck! it fucks everything up :) 
     173!               IF ( ( ABS( d_v_i_thd(ji,jj,jl) ) / MAX( ABS( d_a_i_thd(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 
     174!                  &  .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 
     175!                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp 
     176!                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
     177!               ENDIF 
     178!            END DO ! ji 
     179!         END DO !jj 
     180!      END DO !jl 
    247181      
    248182      at_i(:,:) = 0._wp 
    … …  
    250184         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    251185      END DO 
    252  
    253 !      !IF( ln_nicep ) THEN   
    254 !         WRITE(numout,*) ' limupdate2 - after h correction 2 ' 
    255 !         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
    256 !         WRITE(numout,*) ' at_i : ', at_i(12,44) 
    257 !         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
    258 !         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
    259 !      !ENDIF 
    260186 
    261187      !---------------------------------------------------- 
    … …  
    304230                     zesum = zesum + e_i(ji,jj,jk,jl) 
    305231                  END DO 
    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 
    309232                  ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) 
    310233                  v_i(ji,jj,jl)  = ht_i(ji,jj,jl) * a_i(ji,jj,jl) 
    … …  
    368291            DO ji = 1, jpi 
    369292               ! snow energy of melting 
    370                ze_s = e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), 1.0e-6 )  ! snow energy of melting 
     293               zinda   =  MAX( 0._wp, SIGN( 1._wp, v_s(ji,jj,jl) - epsi10 ) ) 
     294               ze_s = zinda * e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), epsi10 )  ! snow energy of melting 
    371295 
    372296               ! If snow energy of melting smaller then Lf 
    … …  
    401325 
    402326               !switches 
    403                zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
     327               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
    404328               !switch = 1 if a_i > 1e-06 and 0 if not 
    405                zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not 
    406                zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not 
     329               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 
     330               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 
    407331               ! bug fix 25 avril 2007 
    408332               zindb         = zindb*zindic 
    … …  
    444368               !--- 2.7 Correction to ice concentrations  
    445369               !-------------------------------------------- 
    446                !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 ) 
    447370               a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 
    448371 
    … …  
    464387            DO jj = 1, jpj  
    465388               DO ji = 1, jpi 
    466                   zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) )  
     389                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )  
    467390                  e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 
    468391               END DO ! ji 
    … …  
    478401               !--- 2.12 Constrain the thickness of the smallest category above 5 cm 
    479402               !---------------------------------------------------------------------- 
    480                zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
    481                ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi06) 
    482                zh            = MAX( rone , zindb * hiclim  / MAX( ht_i(ji,jj,jl) , epsi20 ) ) 
     403               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
     404               ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi10) 
     405               zh            = MAX( rone , zindb * hiclim  / MAX( ht_i(ji,jj,jl) , epsi10 ) ) 
    483406               ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh 
    484407               ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh 
    … …  
    502425         DO ji = 1, jpi 
    503426            z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )         
    504             zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) )  
     427            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) )  
    505428            DO jl  = 1, jpl 
    506                z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb 
     429               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 
    507430               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 
    508431               ! 
    509                zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
    510                ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda 
     432               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )  
     433               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda 
    511434               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 
    512435            END DO 
    … …  
    542465                     !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) ) 
    543466                     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) 
     467                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 
     468                     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) 
    546469                  END DO ! ji 
    547470               END DO ! jj 
    … …  
    663586            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update2  : old_oa_i    : ') 
    664587            CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update2  : d_oa_i_thd  : ') 
    665             CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update2  : Path. case  : ') 
    666588 
    667589            DO jk = 1, nlay_i 

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

    r4220 r4332  
    6666   PUBLIC   lim_var_salprof1d    ! 
    6767 
    68    REAL(wp) ::   epsi20 = 1.e-20_wp   ! module constants 
    69    REAL(wp) ::   epsi16 = 1.e-16_wp   !    -       - 
    70    REAL(wp) ::   epsi13 = 1.e-13_wp   !    -       - 
    7168   REAL(wp) ::   epsi10 = 1.e-10_wp   !    -       - 
    72    REAL(wp) ::   epsi06 = 1.e-06_wp   !    -       - 
    7369   REAL(wp) ::   zzero = 0.e0        !    -       - 
    7470   REAL(wp) ::   zone  = 1.e0        !    -       - 
    … …  
    117113               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 
    118114               ! 
    119                zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) )  
    120                icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda  ! ice thickness 
     115               zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) )  
     116               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * zinda  ! ice thickness 
    121117            END DO 
    122118         END DO 
    … …  
    138134            DO jj = 1, jpj 
    139135               DO ji = 1, jpi 
    140                   zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi16 ) )  
    141                   zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) )  
     136                  zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi10 ) )  
     137                  zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi10 ) )  
    142138                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content 
    143                   smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi16 ) * zinda   ! ice salinity 
    144                   ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi16 ) * zindb   ! ice age 
     139                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * zinda   ! ice salinity 
     140                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi10 ) * zindb   ! ice age 
    145141               END DO 
    146142            END DO 
    … …  
    209205               DO ji = 1, jpi 
    210206                  !                                                              ! Energy of melting q(S,T) [J.m-3] 
    211                   zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi06 ) * REAL(nlay_i,wp)  
    212                   zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes 
     207                  zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp)  
     208                  zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! zindb = 0 if no ice and 1 if yes 
    213209                  zq_i    = zq_i * unit_fac * zindb                              !convert units 
    214210                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature 
    … …  
    235231               DO ji = 1, jpi 
    236232                  !Energy of melting q(S,T) [J.m-3] 
    237                   zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL(nlay_s,wp) 
    238                   zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes 
     233                  zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp) 
     234                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) )     ! zindb = 0 if no ice and 1 if yes 
    239235                  zq_s  = zq_s * unit_fac * zindb                                    ! convert units 
    240236                  ! 
    … …  
    341337                  zind01 = ( 1._wp - zind0 ) * MAX( 0._wp   , SIGN( 1._wp  , s_i_1 - sm_i(ji,jj,jl) ) )  
    342338                  ! If 2.sm_i GE sss_m then zindbal = 1 
     339                  ! this is to force a constant salinity profile in the Baltic Sea 
    343340                  zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 
    344341                  zalpha(ji,jj,jl) = zind0  + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 ) 
    … …  
    434431            DO jj = 1, jpj 
    435432               DO ji = 1, jpi 
    436                   zinda = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi16 ) )  ) 
    437                   zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi16 ) )  ) 
    438                   zbvi  = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi16 )   & 
     433                  zinda = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) )  ) 
     434                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) 
     435                  zbvi  = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   & 
    439436                     &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp) 
    440                   bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi  / MAX( vt_i(ji,jj) , epsi16 ) 
     437                  bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi  / MAX( vt_i(ji,jj) , epsi10 ) 
    441438               END DO 
    442439            END DO 
    … …  
    498495               zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) )  
    499496               ! if 2.sm_i GE sss_m then zindbal = 1 
     497               ! this is to force a constant salinity profile in the Baltic Sea 
    500498               zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(ii,ij) ) ) 
    501499               ! 

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

    r4045 r4332  
    5252   INTEGER            , DIMENSION(jpnoumax) ::   nc  , nca     ! switch for saving field ( = 1 ) or not ( = 0 ) 
    5353 
    54    REAL(wp)  ::   epsi06 = 1e-6_wp 
     54   REAL(wp)  ::   epsi06 = 1.e-6_wp 
    5555   REAL(wp)  ::   zzero  = 0._wp 
    5656   REAL(wp)  ::   zone   = 1._wp       
    … …  
    192192 
    193193      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
    194       IF( ln_nicep ) THEN 
    195          WRITE(numout,*) 
    196          WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, numit 
    197          WRITE(numout,*) '~~~~~~~ ' 
    198          WRITE(numout,*) ' kindic = ', kindic 
    199       ENDIF 
     194      !IF( ln_nicep ) THEN 
     195      !   WRITE(numout,*) 
     196      !   WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, numit 
     197      !   WRITE(numout,*) '~~~~~~~ ' 
     198      !   WRITE(numout,*) ' kindic = ', kindic 
     199      !ENDIF 
    200200      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
    201201 
    … …  
    564564      CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    565565      CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     566      CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     567      CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     568      CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     569      CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     570      CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     571      CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     572      CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     573      CALL histdef( kid, "iisfxthd", "Salt flux from thermo"   , ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     574      CALL histdef( kid, "iisfxmec", "Salt flux from dynmics"  , ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     575      CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     576 
    566577 
    567578      !CALL histdef( kid, "iice_itd", "Ice concentration by cat", "%"      , jpi, jpj, kh_i, jpl, 1, jpl, -99, 32, "inst(x)", rdt, rdt )  
    … …  
    586597      CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) ) 
    587598 
     599      CALL histwrite( kid, "iicebopr", kt, diag_bot_gr        , jpi*jpj, (/1/) ) 
     600      CALL histwrite( kid, "iicedypr", kt, diag_dyn_gr        , jpi*jpj, (/1/) ) 
     601      CALL histwrite( kid, "iicelapr", kt, diag_lat_gr        , jpi*jpj, (/1/) ) 
     602      CALL histwrite( kid, "iicesipr", kt, diag_sni_gr        , jpi*jpj, (/1/) ) 
     603      CALL histwrite( kid, "iicerepr", kt, diag_res_pr        , jpi*jpj, (/1/) ) 
     604      CALL histwrite( kid, "iicebome", kt, diag_bot_me        , jpi*jpj, (/1/) ) 
     605      CALL histwrite( kid, "iicesume", kt, diag_sur_me        , jpi*jpj, (/1/) ) 
     606      CALL histwrite( kid, "iisfxthd", kt, sfx_thd        , jpi*jpj, (/1/) ) 
     607      CALL histwrite( kid, "iisfxmec", kt, sfx_mec        , jpi*jpj, (/1/) ) 
     608      CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) ) 
     609 
    588610      !CALL histwrite( kid, "iice_itd", kt, a_i  , jpi*jpj*jpl, (/1/)  )   ! area 
    589611      !CALL histwrite( kid, "iice_hid", kt, ht_i , jpi*jpj*jpl, (/1/)  )   ! thickness 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r3651 r4332  
    88   !!            3.3  !  2010-09  (D. Storkey) add ice boundary conditions 
    99   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
     10   !!             -   !  2012-01  (C. Rousset) add ice boundary conditions for lim3 
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_bdy  
    … …  
    4546      REAL, POINTER, DIMENSION(:)     ::  hicif 
    4647      REAL, POINTER, DIMENSION(:)     ::  hsnif 
     48#elif defined key_lim3 
     49      REAL, POINTER, DIMENSION(:,:)   ::  a_i   !: now ice leads fraction climatology 
     50      REAL, POINTER, DIMENSION(:,:)   ::  ht_i  !: Now ice  thickness climatology 
     51      REAL, POINTER, DIMENSION(:,:)   ::  ht_s  !: now snow thickness 
    4752#endif 
    4853   END TYPE OBC_DATA 
    … …  
    7883   REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days 
    7984 
    80 #if defined key_lim2 
    81    INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2              ! Choice of boundary condition for sea ice variables  
    82    INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2_dta          !: = 0 use the initial state as bdy dta ;  
     85#if ( defined key_lim2 || defined key_lim3 ) 
     86   INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim               ! Choice of boundary condition for sea ice variables  
     87   INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim_dta           !: = 0 use the initial state as bdy dta ;  
    8388                                                            !: = 1 read it in a NetCDF file 
    8489#endif 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r3909 r4332  
    1111   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions 
    1212   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
     13   !!             -   !  2012-01  (C. Rousset) add ice boundary conditions for lim3 
    1314   !!---------------------------------------------------------------------- 
    1415#if defined key_bdy 
    … …  
    3132#if defined key_lim2 
    3233   USE ice_2 
     34#elif defined key_lim3 
     35   USE par_ice 
     36   USE ice 
     37   USE limcat_1D          ! redistribute ice input into categories 
    3338#endif 
    3439   USE sbcapr 
    … …  
    4954 
    5055   TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr   ! array of pointers to nbmap 
     56 
     57#if defined key_lim3 
     58   LOGICAL :: ll_bdylim3                  ! determine whether ice input is lim2 (F) or lim3 (T) type 
     59   INTEGER :: jfld_hti, jfld_hts, jfld_ai ! indices of ice thickness, snow thickness and concentration in bf structure 
     60#endif 
    5161 
    5262#  include "domzgr_substitute.h90" 
    … …  
    7787                                                        ! etc. 
    7888      !! 
    79       INTEGER     ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd  ! local indices 
     89      INTEGER     ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl  ! local indices 
    8090      INTEGER,          DIMENSION(jpbgrd) ::   ilen1  
    8191      INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts 
    … …  
    164174 
    165175#if defined key_lim2 
    166             IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN  
     176            IF( nn_ice_lim(ib_bdy) .gt. 0 .and. nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN  
    167177               ilen1(:) = nblen(:) 
    168178               igrd = 1                       ! Everything is at T-points here 
    … …  
    170180                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    171181                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    172                   dta_bdy(ib_bdy)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1)          
     182                  dta_bdy(ib_bdy)%frld (ib) = frld(ii,ij) * tmask(ii,ij,1)          
    173183                  dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1)          
    174184                  dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1)          
    175185               END DO  
     186            ENDIF 
     187#elif defined key_lim3 
     188            IF( nn_ice_lim(ib_bdy) .gt. 0 .and. nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN  
     189               ilen1(:) = nblen(:) 
     190               igrd = 1                       ! Everything is at T-points here 
     191               DO jl = 1, jpl 
     192                  DO ib = 1, ilen1(igrd) 
     193                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     194                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     195                     dta_bdy(ib_bdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1)  
     196                     dta_bdy(ib_bdy)%ht_i(ib,jl) = ht_i(ii,ij,jl) * tmask(ii,ij,1)  
     197                     dta_bdy(ib_bdy)%ht_s(ib,jl) = ht_s(ii,ij,jl) * tmask(ii,ij,1)  
     198                  END DO 
     199               END DO 
    176200            ENDIF 
    177201#endif 
    … …  
    319343                  ENDIF 
    320344               ENDIF 
     345#if defined key_lim3 
     346               IF( .NOT. ll_bdylim3 .AND. nn_ice_lim(ib_bdy) > 0 .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is lim2 type) 
     347                CALL lim_cat_1D ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
     348                                  & dta_bdy(ib_bdy)%ht_i,     dta_bdy(ib_bdy)%ht_s,     dta_bdy(ib_bdy)%a_i     ) 
     349               ENDIF 
     350#endif 
     351 
    321352            ENDIF 
    322353            jstart = jend+1 
    … …  
    366397      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   igrid         ! index for grid type (1,2,3 = T,U,V) 
    367398      INTEGER, POINTER, DIMENSION(:)         ::   nblen, nblenrim  ! short cuts 
     399#if defined key_lim3 
     400      INTEGER, DIMENSION(3) ::   zdimsz   ! number of elements in each of the 4 dimensions (i.e. i,j,t,ice-cat) for an array 
     401      INTEGER               ::   zndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat) 
     402      INTEGER               ::   inum,id1 ! local integer 
     403#endif 
    368404      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures 
    369405      TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   !  
    370406      TYPE(FLD_N) ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read 
    371407#if defined key_lim2 
    372       TYPE(FLD_N) ::   bn_frld, bn_hicif, bn_hsnif      ! 
     408      TYPE(FLD_N) ::   bn_frld, bn_hicif, bn_hsnif       
     409#elif defined key_lim3 
     410      TYPE(FLD_N) ::   bn_a_i, bn_ht_i, bn_ht_s       
    373411#endif 
    374412      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
    375413#if defined key_lim2 
    376414      NAMELIST/nambdy_dta/ bn_frld, bn_hicif, bn_hsnif 
     415#elif defined key_lim3 
     416      NAMELIST/nambdy_dta/ bn_a_i, bn_ht_i, bn_ht_s 
    377417#endif 
    378418      NAMELIST/nambdy_dta/ ln_full_vel 
    … …  
    388428      ! Set nn_dta 
    389429      DO ib_bdy = 1, nb_bdy 
    390          nn_dta(ib_bdy) = MAX(  nn_dyn2d_dta(ib_bdy)       & 
    391                                ,nn_dyn3d_dta(ib_bdy)       & 
    392                                ,nn_tra_dta(ib_bdy)         & 
    393 #if defined key_lim2 
    394                                ,nn_ice_lim2_dta(ib_bdy)    & 
     430         nn_dta(ib_bdy) = MAX( nn_dyn2d_dta(ib_bdy)       & 
     431                              ,nn_dyn3d_dta(ib_bdy)       & 
     432                              ,nn_tra_dta(ib_bdy)         & 
     433#if ( defined key_lim2 || defined key_lim3 ) 
     434                              ,nn_ice_lim_dta(ib_bdy)    & 
    395435#endif 
    396436                              ) 
    … …  
    412452            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
    413453         ENDIF 
    414 #if defined key_lim2 
    415          IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1  ) THEN 
     454#if ( defined key_lim2 || defined key_lim3 ) 
     455         IF( nn_ice_lim(ib_bdy) .gt. 0 .and. nn_ice_lim_dta(ib_bdy) .eq. 1  ) THEN 
    416456            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
    417457         ENDIF 
    … …  
    448488            ln_full_vel = .false. 
    449489            ! ... default values (NB: frequency positive => hours, negative => months) 
    450             !                    !  file       ! frequency !  variable   ! time intep !  clim   ! 'yearly' or ! weights  ! rotation  ! 
    451             !                    !  name       ! hours !   name     !  (T/F)  !  (T/F)  !  'monthly'  ! filename ! pairs     ! 
    452             bn_ssh     = FLD_N(  'bdy_ssh'     ,  24   , 'sossheig' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
    453             bn_u2d     = FLD_N(  'bdy_vel2d_u' ,  24   , 'vobtcrtx' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
    454             bn_v2d     = FLD_N(  'bdy_vel2d_v' ,  24   , 'vobtcrty' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
    455             bn_u3d     = FLD_N(  'bdy_vel3d_u' ,  24   , 'vozocrtx' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
    456             bn_v3d     = FLD_N(  'bdy_vel3d_v' ,  24   , 'vomecrty' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
    457             bn_tem     = FLD_N(  'bdy_tem'     ,  24   , 'votemper' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
    458             bn_sal     = FLD_N(  'bdy_sal'     ,  24   , 'vosaline' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     490            !                 !  file       ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation  ! 
     491            !                 !  name       ! hours     !   name     !  (T/F)     !  (T/F)  !  'monthly'  ! filename ! pairs     ! 
     492            bn_ssh   = FLD_N(  'bdy_ssh'     ,  24      , 'sossheig' , .false.    , .false. ,   'yearly'  , ''       , ''        ) 
     493            bn_u2d   = FLD_N(  'bdy_vel2d_u' ,  24      , 'vobtcrtx' , .false.    , .false. ,   'yearly'  , ''       , ''        ) 
     494            bn_v2d   = FLD_N(  'bdy_vel2d_v' ,  24      , 'vobtcrty' , .false.    , .false. ,   'yearly'  , ''       , ''        ) 
     495            bn_u3d   = FLD_N(  'bdy_vel3d_u' ,  24      , 'vozocrtx' , .false.    , .false. ,   'yearly'  , ''       , ''        ) 
     496            bn_v3d   = FLD_N(  'bdy_vel3d_v' ,  24      , 'vomecrty' , .false.    , .false. ,   'yearly'  , ''       , ''        ) 
     497            bn_tem   = FLD_N(  'bdy_tem'     ,  24      , 'votemper' , .false.    , .false. ,   'yearly'  , ''       , ''        ) 
     498            bn_sal   = FLD_N(  'bdy_sal'     ,  24      , 'vosaline' , .false.    , .false. ,   'yearly'  , ''       , ''        ) 
    459499#if defined key_lim2 
    460             bn_frld    = FLD_N(  'bdy_frld'    ,  24   , 'ildsconc' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
    461             bn_hicif   = FLD_N(  'bdy_hicif'   ,  24   , 'iicethic' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
    462             bn_hsnif   = FLD_N(  'bdy_hsnif'   ,  24   , 'isnothic' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     500            bn_frld  = FLD_N(  'bdy_frld'    ,  24      , 'ildsconc' , .false.    , .false. ,   'yearly'  , ''       , ''        ) 
     501            bn_hicif = FLD_N(  'bdy_hicif'   ,  24      , 'iicethic' , .false.    , .false. ,   'yearly'  , ''       , ''        ) 
     502            bn_hsnif = FLD_N(  'bdy_hsnif'   ,  24      , 'isnothic' , .false.    , .false. ,   'yearly'  , ''       , ''        ) 
     503#elif defined key_lim3 
     504            bn_a_i  = FLD_N(  'bdy_a_i'     ,  24   , 'ileadfra' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     505            bn_ht_i = FLD_N(  'bdy_ht_i'    ,  24   , 'iicethic' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
     506            bn_ht_s = FLD_N(  'bdy_ht_s'    ,  24   , 'isnowthi' , .false. , .false. ,   'yearly'  , ''       , ''        ) 
    463507#endif 
    464508 
    … …  
    545589#if defined key_lim2 
    546590            ! sea ice 
    547             IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN 
     591            IF( nn_ice_lim(ib_bdy) .gt. 0 .and. nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 
    548592 
    549593               jfld = jfld + 1 
    … …  
    567611               ilen1(jfld) = nblen(igrid(jfld)) 
    568612               ilen3(jfld) = 1 
     613 
     614            ENDIF 
     615#elif defined key_lim3 
     616            ! sea ice 
     617            IF( nn_ice_lim(ib_bdy) .gt. 0 .and. nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 
     618 
     619               ! Test for types of ice input (lim2 or lim3)  
     620               CALL iom_open ( bn_a_i%clname, inum ) 
     621               id1 = iom_varid ( inum, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 
     622               CALL iom_close ( inum ) 
     623               !CALL fld_clopn ( bn_a_i, nyear, nmonth, nday, ldstop=.TRUE. ) 
     624               !CALL iom_open ( bn_a_i%clname, inum ) 
     625               !id1 = iom_varid ( bn_a_i%num, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 
     626                IF ( zndims == 4 ) THEN 
     627                 ll_bdylim3 = .TRUE.   ! lim3 input 
     628               ELSE 
     629                 ll_bdylim3 = .FALSE.  ! lim2 input       
     630               ENDIF 
     631               ! End test 
     632 
     633               jfld = jfld + 1 
     634               blf_i(jfld) = bn_a_i 
     635               ibdy(jfld) = ib_bdy 
     636               igrid(jfld) = 1 
     637               ilen1(jfld) = nblen(igrid(jfld)) 
     638               IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 
     639 
     640               jfld = jfld + 1 
     641               blf_i(jfld) = bn_ht_i 
     642               ibdy(jfld) = ib_bdy 
     643               igrid(jfld) = 1 
     644               ilen1(jfld) = nblen(igrid(jfld)) 
     645               IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 
     646 
     647               jfld = jfld + 1 
     648               blf_i(jfld) = bn_ht_s 
     649               ibdy(jfld) = ib_bdy 
     650               igrid(jfld) = 1 
     651               ilen1(jfld) = nblen(igrid(jfld)) 
     652               IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 
    569653 
    570654            ENDIF 
    … …  
    662746 
    663747#if defined key_lim2 
    664          IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN 
    665             IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
     748         IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 
     749            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 
    666750               ilen0(1:3) = nblen(1:3) 
    667751               ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) ) 
    … …  
    675759               jfld = jfld + 1 
    676760               dta_bdy(ib_bdy)%hsnif => bf(jfld)%fnow(:,1,1) 
     761            ENDIF 
     762         ENDIF 
     763#elif defined key_lim3 
     764         IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 
     765            ilen0(1:3) = nblen(1:3) 
     766            IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 
     767               ALLOCATE( dta_bdy(ib_bdy)%a_i (ilen0(1),jpl) ) 
     768               ALLOCATE( dta_bdy(ib_bdy)%ht_i(ilen0(1),jpl) ) 
     769               ALLOCATE( dta_bdy(ib_bdy)%ht_s(ilen0(1),jpl) ) 
     770            ELSE 
     771               IF ( ll_bdylim3 ) THEN ! case input is lim3 type 
     772                  jfld = jfld + 1 
     773                  dta_bdy(ib_bdy)%a_i  => bf(jfld)%fnow(:,1,:) 
     774                  jfld = jfld + 1 
     775                  dta_bdy(ib_bdy)%ht_i => bf(jfld)%fnow(:,1,:) 
     776                  jfld = jfld + 1 
     777                  dta_bdy(ib_bdy)%ht_s => bf(jfld)%fnow(:,1,:) 
     778               ELSE ! case input is lim2 type 
     779                  jfld_ai  = jfld + 1 
     780                  jfld_hti = jfld + 2 
     781                  jfld_hts = jfld + 3 
     782                  jfld     = jfld + 3 
     783                  ALLOCATE( dta_bdy(ib_bdy)%a_i (ilen0(1),jpl) ) 
     784                  ALLOCATE( dta_bdy(ib_bdy)%ht_i(ilen0(1),jpl) ) 
     785                  ALLOCATE( dta_bdy(ib_bdy)%ht_s(ilen0(1),jpl) ) 
     786                  dta_bdy(ib_bdy)%a_i (:,:) = 0.0 
     787                  dta_bdy(ib_bdy)%ht_i(:,:) = 0.0 
     788                  dta_bdy(ib_bdy)%ht_s(:,:) = 0.0 
     789               ENDIF 
     790 
    677791            ENDIF 
    678792         ENDIF 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r3909 r4332  
    9696         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         &   
    9797         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp,              & 
    98 #if defined key_lim2 
    99          &             nn_ice_lim2, nn_ice_lim2_dta,                       & 
     98#if ( defined key_lim2 || defined key_lim3 ) 
     99         &             nn_ice_lim, nn_ice_lim_dta,                       & 
    100100#endif 
    101101         &             ln_vol, nn_volctl, nn_rimwidth 
    … …  
    137137      ln_dyn3d_dmp(:)   = .false. 
    138138      rn_time_dmp(:)    = 1. 
    139 #if defined key_lim2 
    140       nn_ice_lim2(:)    = 0 
    141       nn_ice_lim2_dta(:)= -1  ! uninitialised flag 
     139#if ( defined key_lim2 || defined key_lim3 ) 
     140      nn_ice_lim(:)    = 0 
     141      nn_ice_lim_dta(:)= -1  ! uninitialised flag 
    142142#endif 
    143143      ln_vol            = .false. 
    … …  
    161161 
    162162      DO ib_bdy = 1,nb_bdy 
     163        icount = 0 ! flag to set max rimwidth to 1 if no relaxation 
    163164        IF(lwp) WRITE(numout,*) ' '  
    164165        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'  
    … …  
    175176          CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    176177          CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     178          icount = icount + 1 
    177179          CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
    178180          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 
    … …  
    196198          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    197199          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     200          icount = icount + 1 
    198201          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    199202          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
    … …  
    215218              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    216219           ELSE 
     220              icount = icount + 1 
    217221              IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone' 
    218222              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
    … …  
    228232          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    229233          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     234          icount = icount + 1 
    230235          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    231236          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions' 
    … …  
    248253              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    249254           ELSE 
     255              icount = icount + 1 
    250256              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
    251257              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
    … …  
    257263        IF(lwp) WRITE(numout,*) 
    258264 
    259 #if defined key_lim2 
     265#if ( defined key_lim2 || defined key_lim3 ) 
    260266        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  ' 
    261         SELECT CASE( nn_ice_lim2(ib_bdy) )                   
     267        SELECT CASE( nn_ice_lim(ib_bdy) )                   
    262268          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    263269          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     270          icount = icount + 1 
    264271          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
    265272        END SELECT 
    266         IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN  
    267            SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   !  
     273        IF( nn_ice_lim(ib_bdy) .gt. 0 ) THEN  
     274           SELECT CASE( nn_ice_lim_dta(ib_bdy) )                   !  
    268275              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    269276              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file' 
    270               CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' ) 
     277              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim_dta must be 0 or 1' ) 
    271278           END SELECT 
    272279        ENDIF 
    273280        IF(lwp) WRITE(numout,*) 
    274281#endif 
    275  
    276         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
    277         IF(lwp) WRITE(numout,*) 
     282        IF ( icount>0 ) THEN 
     283           IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
     284           IF(lwp) WRITE(numout,*) 
     285        ELSE 
     286           nn_rimwidth(ib_bdy) = 1 ! no relaxation 
     287        ENDIF 
    278288 
    279289      ENDDO 
    … …  
    387397            DO igrd = 1, jpbgrd 
    388398               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
    389                nblendta(igrd,ib_bdy) = kdimsz(1) 
    390                jpbdtau = MAX(jpbdtau, kdimsz(1)) 
     399               !clem nblendta(igrd,ib_bdy) = kdimsz(1) 
     400               !clem jpbdtau = MAX(jpbdtau, kdimsz(1)) 
     401               nblendta(igrd,ib_bdy) = MAXVAL(kdimsz) 
     402               jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz)) 
    391403            ENDDO 
    392404            CALL iom_close( inum ) 
    393  
    394405         ENDIF  
    395406 
    … …  
    398409      IF (nb_bdy>0) THEN 
    399410         jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 
    400  
    401411         ! Allocate arrays 
    402412         !--------------- 
    … …  
    446456         ENDIF  
    447457 
    448       ENDDO       
     458      ENDDO      
    449459     
    450460      ! 2. Now fill indices corresponding to straight open boundary arrays: 
    … …  
    752762               ! check if point is in local domain 
    753763               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   & 
    754                   & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN 
     764                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   & 
     765                  & nbrdta(ib,igrd,ib_bdy) <= nn_rimwidth(ib_bdy)     ) THEN       
    755766                  ! 
    756767                  icount = icount  + 1 
    … …  
    765776         ! Allocate index arrays for this boundary set 
    766777         !-------------------------------------------- 
    767          ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:)) 
     778 
     779         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(1:jpbgrd)) 
     780         ilen1 = MAX(1,ilen1) 
    768781         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) ) 
    769782         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 
    … …  
    773786         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 
    774787         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) ) 
    775          ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) ) 
     788         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )       
    776789 
    777790         ! Dispatch mapping indices and discrete distances on each processor 
    … …  
    10081021      ! bdytmask = 1  on the computational domain AND on open boundaries 
    10091022      !          = 0  elsewhere    
    1010   
     1023      bdytmask(:,:) = 1.e0 
     1024      bdyumask(:,:) = 1.e0 
     1025      bdyvmask(:,:) = 1.e0 
     1026 
    10111027      IF( ln_mask_file ) THEN 
    10121028         CALL iom_open( cn_mask_file, inum ) 
    … …  
    10531069       
    10541070      bdytmask(:,:) = tmask(:,:,1) 
    1055       IF( .not. ln_mask_file ) THEN 
    1056          ! If .not. ln_mask_file then we need to derive mask on U and V grid  
    1057          ! from mask on T grid here. 
    1058          bdyumask(:,:) = 0.e0 
    1059          bdyvmask(:,:) = 0.e0 
    1060          DO ij=1, jpjm1 
    1061             DO ii=1, jpim1 
    1062                bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 
    1063                bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)   
    1064             END DO 
    1065          END DO 
    1066          CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond. 
    1067       ENDIF 
    10681071 
    10691072      ! bdy masks and bmask are now set to zero on boundary points: 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r4220 r4332  
    237237      ENDIF 
    238238      ! 
    239       IF( lk_lim3 ) THEN  
     239      IF( lk_lim3 ) THEN 
    240240         CALL iom_get( numror, jpdom_autoglo, 'iatte' , iatte ) ! clem modif 
    241241         CALL iom_get( numror, jpdom_autoglo, 'oatte' , oatte ) ! clem modif 

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

    r4040 r4332  
    724724        Cd_n10(:,:) =   cdn_wave 
    725725      ELSE 
    726         Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a) 
     726        Cd_n10  = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a) 
    727727      ENDIF 
    728728      sqrt_Cd_n10 = sqrt(Cd_n10) 
    729       Ce_n10  = 1E-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b) 
    730       Ch_n10  = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1-stab)) !   L & Y eq. (6c), (6d) 
     729      Ce_n10  = 1.e-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b) 
     730      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) !   L & Y eq. (6c), (6d) 
    731731      !! 
    732732      !! Initializing transfert coefficients with their first guess neutral equivalents : 
    … …  
    755755 
    756756           !! Updating the neutral 10m transfer coefficients : 
    757            Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a) 
     757           Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a) 
    758758           sqrt_Cd_n10 = sqrt(Cd_n10) 
    759            Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b) 
     759           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b) 
    760760           stab    = 0.5 + sign(0.5,zeta) 
    761            Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab))           !  L & Y eq. (6c), (6d) 
     761           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab))           !  L & Y eq. (6c), (6d) 
    762762 
    763763           !! Shifting the neutral  10m transfer coefficients to ( zu , zeta ) : 
    … …  
    858858        Cd_n10(:,:) =   cdn_wave 
    859859      ELSE 
    860         Cd_n10  = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 )  
     860        Cd_n10  = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 )  
    861861      ENDIF 
    862862      sqrt_Cd_n10 = sqrt(Cd_n10) 
    863       Ce_n10  = 1E-3*( 34.6 * sqrt_Cd_n10 ) 
    864       Ch_n10  = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1 - stab)) 
     863      Ce_n10  = 1.e-3*( 34.6 * sqrt_Cd_n10 ) 
     864      Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab)) 
    865865 
    866866      !! Initializing transf. coeff. with their first guess neutral equivalents : 
    … …  
    904904         ELSE 
    905905           !! Updating the neutral 10m transfer coefficients : 
    906            Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a) 
     906           Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)    ! L & Y eq. (6a) 
    907907           sqrt_Cd_n10 = sqrt(Cd_n10) 
    908            Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b) 
     908           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                 ! L & Y eq. (6b) 
    909909           stab    = 0.5 + sign(0.5,zeta_u) 
    910            Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d) 
     910           Ch_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1.-stab)) ! L & Y eq. (6c-6d) 
    911911           !! 
    912912           !! 

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

    r4220 r4332  
    2727   USE iceini          ! LIM-3: ice initialisation 
    2828   USE dom_ice         ! LIM-3: ice domain 
     29   USE domvvl          ! domain: variable volume level(fse3t_m) 
    2930 
    3031   USE sbc_oce         ! Surface boundary condition: ocean fields 
    … …  
    5859   USE in_out_manager  ! I/O manager 
    5960   USE prtctl          ! Print control 
     61   USE lib_fortran     !  
    6062 
    6163#if defined key_bdy  
    … …  
    167169         ! 
    168170         IF( ln_nicep ) THEN      ! control print at a given point 
    169             jiindx = 6   ;   jjindx = 47 
    170             WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
     171            jiindx = 177   ;   jjindx = 112 
     172            IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    171173         ENDIF 
    172174      ENDIF 
    … …  
    303305         fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp 
    304306         fhmec  (:,:) = 0._wp   ;    
    305          fmmec  (:,:) = 0._wp      
     307         fmmec  (:,:) = 0._wp 
     308         fmmflx (:,:) = 0._wp      
    306309         focea2D(:,:) = 0._wp 
    307310         fsup2D (:,:) = 0._wp 
    … …  
    328331                          CALL lim_rst_opn( kt )     ! Open Ice restart file 
    329332         ! 
    330          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
    331          ! 
    332          IF( .NOT. lk_c1d ) THEN                     ! Ice dynamics & transport (except in 1D case) 
     333         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
     334         ! ---------------------------------------------- 
     335         ! ice dynamics and transport (except in 1D case) 
     336         ! ---------------------------------------------- 
     337         IF( .NOT. lk_c1d ) THEN 
    333338                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics ) 
    334339                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion ) 
    335340                          CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
    336          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
     341         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
    337342                          CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting) 
    338343                          CALL lim_var_agg( 1 )  
     344#if defined key_bdy 
     345                          ! bdy ice thermo  
     346                          CALL lim_var_glo2eqv            ! equivalent variables 
     347                          CALL bdy_ice_lim( kt ) 
     348                          CALL lim_itd_me_zapsmall 
     349                          CALL lim_var_agg(1) 
     350         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print 
     351#endif 
    339352                          CALL lim_update1 
    340353         ENDIF 
    … …  
    349362                          old_oa_i(:,:,:)  = oa_i(:,:,:) 
    350363                          old_smv_i(:,:,:) = smv_i (:,:,:) 
    351          !                                           ! Ice thermodynamics  
     364  
     365         ! ---------------------------------------------- 
     366         ! ice thermodynamic 
     367         ! ---------------------------------------------- 
    352368                          CALL lim_var_glo2eqv            ! equivalent variables 
    353369                          CALL lim_var_agg(1)             ! aggregate ice categories 
     370                          ! previous lead fraction and ice volume for flux calculations 
     371                          pfrld(:,:)   = 1._wp - at_i(:,:) 
     372                          phicif(:,:)  = vt_i(:,:) 
     373                          ! 
    354374                          CALL lim_var_bv                 ! bulk brine volume (diag) 
    355375                          CALL lim_thd( kt )              ! Ice thermodynamics  
    … …  
    357377                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 
    358378                          CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin) 
    359          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
     379         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
    360380                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  ! 
    361          ! 
    362          !                                           ! Global variables update 
    363381                          CALL lim_var_agg( 1 )           ! requested by limupdate 
    364                           CALL lim_update2                 ! Global variables update 
    365 #if defined key_bdy 
    366                           CALL lim_var_glo2eqv            ! 
    367                           CALL bdy_ice_lim( kt )          ! clem modif: bdy ice 
    368 #endif 
     382                          CALL lim_update2                ! Global variables update 
     383 
    369384                          CALL lim_var_glo2eqv            ! equivalent variables (outputs) 
    370385                          CALL lim_var_agg(2)             ! aggregate ice thickness categories 
    371          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
     386         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
    372387         ! 
    373388                          CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes 
    374389         ! 
    375          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
     390         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
    376391         ! 
    377392         !                                           ! Diagnostics and outputs  
    … …  
    386401                          CALL lim_var_glo2eqv            ! ??? 
    387402         ! 
    388          IF( ln_nicep )   CALL lim_ctl               ! alerts in case of model crash 
     403         IF( ln_nicep )   CALL lim_ctl( kt )              ! alerts in case of model crash 
    389404         ! 
    390405      ENDIF                                    ! End sea-ice time step only 
    … …  
    413428 
    414429 
    415    SUBROUTINE lim_ctl 
     430   SUBROUTINE lim_ctl( kt ) 
    416431      !!----------------------------------------------------------------------- 
    417432      !!                   ***  ROUTINE lim_ctl ***  
    … …  
    419434      !! ** Purpose :   Alerts in case of model crash 
    420435      !!------------------------------------------------------------------- 
     436      INTEGER, INTENT(in) ::   kt      ! ocean time step 
    421437      INTEGER  ::   ji, jj, jk,  jl   ! dummy loop indices 
    422438      INTEGER  ::   inb_altests       ! number of alert tests (max 20) 
    … …  
    438454            DO ji = 1, jpi 
    439455               IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN 
    440                   WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    441                   WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
    442                   WRITE(numout,*) ' Point - category', ji, jj, jl 
    443                   WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl) 
    444                   WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl) 
    445                   WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl) 
    446                   WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl) 
     456                  !WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
     457                  !WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
     458                  !WRITE(numout,*) ' Point - category', ji, jj, jl 
     459                  !WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl) 
     460                  !WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl) 
     461                  !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl) 
     462                  !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl) 
    447463                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    448464               ENDIF 
    … …  
    458474         DO ji = 1, jpi 
    459475            IF(   ht_i(ji,jj,jl)  >  50._wp   ) THEN 
    460                CALL lim_prt_state( ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
     476               !CALL lim_prt_state( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
    461477               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    462478            ENDIF 
    … …  
    469485      DO jj = 1, jpj 
    470486         DO ji = 1, jpi 
    471             IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 0.5  .AND.  & 
     487            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  & 
    472488               &  at_i(ji,jj) > 0._wp   ) THEN 
    473                CALL lim_prt_state( ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
    474                WRITE(numout,*) ' ice strength             : ', strength(ji,jj) 
    475                WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj)  
    476                WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj) 
    477                WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)  
    478                WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj) 
    479                WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj) 
    480                WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj) 
    481                WRITE(numout,*) ' sst                      : ', sst_m(ji,jj) 
    482                WRITE(numout,*) ' sss                      : ', sss_m(ji,jj) 
    483                WRITE(numout,*)  
     489               !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
     490               !WRITE(numout,*) ' ice strength             : ', strength(ji,jj) 
     491               !WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj)  
     492               !WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj) 
     493               !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)  
     494               !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj) 
     495               !WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj) 
     496               !WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj) 
     497               !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj) 
     498               !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj) 
     499               !WRITE(numout,*)  
    484500               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    485501            ENDIF 
    … …  
    493509         DO ji = 1, jpi 
    494510            IF(   tms(ji,jj) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN  
    495                CALL lim_prt_state( ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
    496                WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)  
    497                WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    498                WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    499                WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj) 
    500                WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj) 
    501                WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1) 
    502                WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj) 
    503                WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj) 
     511               !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
     512               !WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)  
     513               !WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
     514               !WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
     515               !WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj) 
     516               !WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj) 
     517               !WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1) 
     518               !WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj) 
     519               !WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj) 
    504520               ! 
    505521               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    … …  
    515531         DO jj = 1, jpj 
    516532            DO ji = 1, jpi 
    517 !!gm  test twice sm_i ...  ????  bug? 
    518                IF( ( ( ABS( sm_i(ji,jj,jl) ) < 0.5 )   .OR.  & 
    519                      ( ABS( sm_i(ji,jj,jl) ) < 0.5 ) ) .AND. & 
    520                              ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
    521 !                 CALL lim_prt_state(ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
     533               IF( sm_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
     534!                 CALL lim_prt_state(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
    522535!                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    523536!                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    … …  
    540553                      ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 
    541554                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
    542                   CALL lim_prt_state( ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
     555                  !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
    543556                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    544557               ENDIF 
    … …  
    552565      DO jj = 1, jpj 
    553566         DO ji = 1, jpi 
    554             IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN 
    555                CALL lim_prt_state( ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    556                DO jl = 1, jpl 
    557                   WRITE(numout,*) ' Category no: ', jl 
    558                   WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)    
    559                   WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    560                   WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)    
    561                   WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    562                   WRITE(numout,*) ' ' 
    563                END DO 
     567            IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
     568               !CALL lim_prt_state( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
     569               !DO jl = 1, jpl 
     570                  !WRITE(numout,*) ' Category no: ', jl 
     571                  !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)    
     572                  !WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
     573                  !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)    
     574                  !WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
     575                  !WRITE(numout,*) ' ' 
     576               !END DO 
    564577               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    565578            ENDIF 
    … …  
    572585      DO jj = 1, jpj 
    573586         DO ji = 1, jpi 
    574             IF(   ABS( qns(ji,jj) ) > 1500._wp  .AND.  at_i(ji,jj) > 0._wp  ) THEN 
     587            IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
    575588               ! 
    576                WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
    577                WRITE(numout,*) ' ji, jj    : ', ji, jj 
    578                WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    579                WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
    580                WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
    581                WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) 
    582                WRITE(numout,*) ' qldif     : ', qldif(ji,jj) 
    583                WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice 
    584                WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice 
    585                WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj) 
    586                WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj) 
    587                WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice 
    588                WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice 
    589                WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)  
    590                WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj)  
    591                WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj)  
    592                WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)  
    593                WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj)  
     589               !WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
     590               !WRITE(numout,*) ' ji, jj    : ', ji, jj 
     591               !WRITE(numout,*) ' qns       : ', qns(ji,jj) 
     592               !WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
     593               !WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
     594               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) 
     595               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) 
     596               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice 
     597               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice 
     598               !WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj) 
     599               !WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj) 
     600               !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice 
     601               !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice 
     602               !WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)  
     603               !WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj)  
     604               !WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj)  
     605               !WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)  
     606               !WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj)  
    594607               ! 
    595                CALL lim_prt_state( ji, jj, 2, '   ') 
     608               !CALL lim_prt_state( kt, ji, jj, 2, '   ') 
    596609               inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    597610               ! 
    … …  
    610623               DO ji = 1, jpi 
    611624                  ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt 
    612                   IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-6   & 
     625                  IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
    613626                     &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
    614                      WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
    615                      WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 
    616                      WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl) 
    617                      WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl) 
    618                      WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl) 
    619                      WRITE(numout,*) ' ztmelts : ', ztmelts 
     627                     !WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
     628                     !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl 
     629                     !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl) 
     630                     !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl) 
     631                     !WRITE(numout,*) ' s_i : ', s_i(ji,jj,jk,jl) 
     632                     !WRITE(numout,*) ' ztmelts : ', ztmelts 
    620633                     inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    621634                  ENDIF 
    … …  
    625638      END DO 
    626639 
    627       ialert_id = 1                                 ! reference number of this alert 
    628       cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert 
    629       WRITE(numout,*) 
    630       WRITE(numout,*) ' All alerts at the end of ice model ' 
    631       DO ialert_id = 1, inb_altests 
    632          WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! ' 
    633       END DO 
    634       ! 
     640      ! sum of the alerts on all processors 
     641      IF( lk_mpp ) THEN 
     642         DO ialert_id = 1, inb_altests 
     643            CALL mpp_sum(inb_alp(ialert_id)) 
     644         END DO 
     645      ENDIF 
     646 
     647      ! print alerts 
     648      IF( lwp ) THEN 
     649         ialert_id = 1                                 ! reference number of this alert 
     650         cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert 
     651         WRITE(numout,*) ' time step ',kt 
     652         WRITE(numout,*) ' All alerts at the end of ice model ' 
     653         DO ialert_id = 1, inb_altests 
     654            WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! ' 
     655         END DO 
     656      ENDIF 
     657     ! 
    635658   END SUBROUTINE lim_ctl 
    636659  
    637660    
    638    SUBROUTINE lim_prt_state( ki, kj, kn, cd1 ) 
     661   SUBROUTINE lim_prt_state( kt, ki, kj, kn, cd1 ) 
    639662      !!----------------------------------------------------------------------- 
    640663      !!                   ***  ROUTINE lim_prt_state ***  
    … …  
    650673      !!                n : number of the option 
    651674      !!------------------------------------------------------------------- 
     675      INTEGER         , INTENT(in) ::   kt      ! ocean time step 
    652676      INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices 
    653677      CHARACTER(len=*), INTENT(in) ::   cd1           ! 
    654678      !! 
    655       INTEGER :: jl 
     679      INTEGER :: jl, ji, jj 
    656680      !!------------------------------------------------------------------- 
    657681 
    658       WRITE(numout,*) cd1             ! print title 
    659  
    660       !---------------- 
    661       !  Simple state 
    662       !---------------- 
    663  
    664       IF ( kn == 1 .OR. kn == -1 ) THEN 
    665          WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj 
    666          WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    667          WRITE(numout,*) ' Simple state ' 
    668          WRITE(numout,*) ' masks s,u,v   : ', tms(ki,kj), tmu(ki,kj), tmv(ki,kj) 
    669          WRITE(numout,*) ' lat - long    : ', gphit(ki,kj), glamt(ki,kj) 
    670          WRITE(numout,*) ' Time step     : ', numit 
    671          WRITE(numout,*) ' - Ice drift   ' 
    672          WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    673          WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ki-1,kj) 
    674          WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ki,kj) 
    675          WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ki,kj-1) 
    676          WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ki,kj) 
    677          WRITE(numout,*) ' strength      : ', strength(ki,kj) 
    678          WRITE(numout,*) 
    679          WRITE(numout,*) ' - Cell values ' 
    680          WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    681          WRITE(numout,*) ' cell area     : ', area(ki,kj) 
    682          WRITE(numout,*) ' at_i          : ', at_i(ki,kj)        
    683          WRITE(numout,*) ' vt_i          : ', vt_i(ki,kj)        
    684          WRITE(numout,*) ' vt_s          : ', vt_s(ki,kj)        
    685          DO jl = 1, jpl 
    686             WRITE(numout,*) ' - Category (', jl,')' 
    687             WRITE(numout,*) ' a_i           : ', a_i(ki,kj,jl) 
    688             WRITE(numout,*) ' ht_i          : ', ht_i(ki,kj,jl) 
    689             WRITE(numout,*) ' ht_s          : ', ht_s(ki,kj,jl) 
    690             WRITE(numout,*) ' v_i           : ', v_i(ki,kj,jl) 
    691             WRITE(numout,*) ' v_s           : ', v_s(ki,kj,jl) 
    692             WRITE(numout,*) ' e_s           : ', e_s(ki,kj,1,jl)/1.0e9 
    693             WRITE(numout,*) ' e_i           : ', e_i(ki,kj,1:nlay_i,jl)/1.0e9 
    694             WRITE(numout,*) ' t_su          : ', t_su(ki,kj,jl) 
    695             WRITE(numout,*) ' t_snow        : ', t_s(ki,kj,1,jl) 
    696             WRITE(numout,*) ' t_i           : ', t_i(ki,kj,1:nlay_i,jl) 
    697             WRITE(numout,*) ' sm_i          : ', sm_i(ki,kj,jl) 
    698             WRITE(numout,*) ' smv_i         : ', smv_i(ki,kj,jl) 
    699             WRITE(numout,*) 
    700             WRITE(numout,*) ' Pathological case : ', patho_case(ki,kj,jl) 
    701          END DO 
    702       ENDIF 
    703       IF( kn == -1 ) THEN 
    704          WRITE(numout,*) ' Mechanical Check ************** ' 
    705          WRITE(numout,*) ' Check what means ice divergence ' 
    706          WRITE(numout,*) ' Total ice concentration ', at_i (ki,kj) 
    707          WRITE(numout,*) ' Total lead fraction     ', ato_i(ki,kj) 
    708          WRITE(numout,*) ' Sum of both             ', ato_i(ki,kj) + at_i(ki,kj) 
    709          WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ki,kj) + at_i(ki,kj) - 1.00 
    710       ENDIF 
    711  
    712  
    713       !-------------------- 
    714       !  Exhaustive state 
    715       !-------------------- 
    716  
    717       IF ( kn .EQ. 2 ) THEN 
    718          WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj 
    719          WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    720          WRITE(numout,*) ' Exhaustive state ' 
    721          WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj) 
    722          WRITE(numout,*) ' Time step ', numit 
    723          WRITE(numout,*)  
    724          WRITE(numout,*) ' - Cell values ' 
    725          WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    726          WRITE(numout,*) ' cell area     : ', area(ki,kj) 
    727          WRITE(numout,*) ' at_i          : ', at_i(ki,kj)        
    728          WRITE(numout,*) ' vt_i          : ', vt_i(ki,kj)        
    729          WRITE(numout,*) ' vt_s          : ', vt_s(ki,kj)        
    730          WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ki-1,kj) 
    731          WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ki,kj) 
    732          WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ki,kj-1) 
    733          WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ki,kj) 
    734          WRITE(numout,*) ' strength      : ', strength(ki,kj) 
    735          WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ki,kj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ki,kj) 
    736          WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ki,kj)  , ' old_v_ice     : ', old_v_ice(ki,kj)   
    737          WRITE(numout,*) 
    738  
    739          DO jl = 1, jpl 
    740               WRITE(numout,*) ' - Category (',jl,')' 
    741               WRITE(numout,*) '   ~~~~~~~~         '  
    742               WRITE(numout,*) ' ht_i       : ', ht_i(ki,kj,jl)             , ' ht_s       : ', ht_s(ki,kj,jl) 
    743               WRITE(numout,*) ' t_i        : ', t_i(ki,kj,1:nlay_i,jl) 
    744               WRITE(numout,*) ' t_su       : ', t_su(ki,kj,jl)             , ' t_s        : ', t_s(ki,kj,1,jl) 
    745               WRITE(numout,*) ' sm_i       : ', sm_i(ki,kj,jl)             , ' o_i        : ', o_i(ki,kj,jl) 
    746               WRITE(numout,*) ' a_i        : ', a_i(ki,kj,jl)              , ' old_a_i    : ', old_a_i(ki,kj,jl)    
    747               WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ki,kj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ki,kj,jl)  
    748               WRITE(numout,*) ' v_i        : ', v_i(ki,kj,jl)              , ' old_v_i    : ', old_v_i(ki,kj,jl)    
    749               WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ki,kj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ki,kj,jl)  
    750               WRITE(numout,*) ' v_s        : ', v_s(ki,kj,jl)              , ' old_v_s    : ', old_v_s(ki,kj,jl)   
    751               WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ki,kj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ki,kj,jl) 
    752               WRITE(numout,*) ' e_i1       : ', e_i(ki,kj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ki,kj,1,jl)/1.0e9  
    753               WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ki,kj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ki,kj,1,jl)/1.0e9 
    754               WRITE(numout,*) ' e_i2       : ', e_i(ki,kj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ki,kj,2,jl)/1.0e9   
    755               WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ki,kj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ki,kj,2,jl)/1.0e9 
    756               WRITE(numout,*) ' e_snow     : ', e_s(ki,kj,1,jl)            , ' old_e_snow : ', old_e_s(ki,kj,1,jl)  
    757               WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ki,kj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ki,kj,1,jl) 
    758               WRITE(numout,*) ' smv_i      : ', smv_i(ki,kj,jl)            , ' old_smv_i  : ', old_smv_i(ki,kj,jl)    
    759               WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ki,kj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ki,kj,jl)  
    760               WRITE(numout,*) ' oa_i       : ', oa_i(ki,kj,jl)             , ' old_oa_i   : ', old_oa_i(ki,kj,jl) 
    761               WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ki,kj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ki,kj,jl) 
    762               WRITE(numout,*) ' Path. case : ', patho_case(ki,kj,jl) 
    763         END DO !jl 
    764  
    765         WRITE(numout,*) 
    766         WRITE(numout,*) ' - Heat / FW fluxes ' 
    767         WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    768 !       WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ki,kj) 
    769 !       WRITE(numout,*) ' sfx        : ', sfx      (ki,kj) 
    770 !       WRITE(numout,*) ' sfx_res  : ', sfx_res(ki,kj) 
    771         WRITE(numout,*) ' fmmec      : ', fmmec    (ki,kj) 
    772         WRITE(numout,*) ' fhmec      : ', fhmec    (ki,kj) 
    773         WRITE(numout,*) ' fhbri      : ', fhbri    (ki,kj) 
    774         WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ki,kj) 
    775         WRITE(numout,*)  
    776         WRITE(numout,*) ' sst        : ', sst_m(ki,kj)   
    777         WRITE(numout,*) ' sss        : ', sss_m(ki,kj)   
    778         WRITE(numout,*)  
    779         WRITE(numout,*) ' - Stresses ' 
    780         WRITE(numout,*) '   ~~~~~~~~ ' 
    781         WRITE(numout,*) ' utau_ice   : ', utau_ice(ki,kj)  
    782         WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ki,kj) 
    783         WRITE(numout,*) ' utau       : ', utau    (ki,kj)  
    784         WRITE(numout,*) ' vtau       : ', vtau    (ki,kj) 
    785         WRITE(numout,*) ' oc. vel. u : ', u_oce   (ki,kj) 
    786         WRITE(numout,*) ' oc. vel. v : ', v_oce   (ki,kj) 
    787      ENDIF 
    788  
    789      !--------------------- 
    790      ! Salt / heat fluxes 
    791      !--------------------- 
    792  
    793      IF ( kn .EQ. 3 ) THEN 
    794         WRITE(numout,*) ' lim_prt_state - Point : ',ki,kj 
    795         WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    796         WRITE(numout,*) ' - Salt / Heat Fluxes ' 
    797         WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    798         WRITE(numout,*) ' lat - long ', gphit(ki,kj), glamt(ki,kj) 
    799         WRITE(numout,*) ' Time step ', numit 
    800         WRITE(numout,*) 
    801         WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
    802         WRITE(numout,*) ' qsr       : ', qsr(ki,kj) 
    803         WRITE(numout,*) ' qns       : ', qns(ki,kj) 
    804         WRITE(numout,*) 
    805         WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
    806         WRITE(numout,*) ' emp       : ', emp    (ki,kj) 
    807         WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ki,kj) 
    808         WRITE(numout,*) ' sfx       : ', sfx    (ki,kj) 
    809         WRITE(numout,*) ' sfx_res   : ', sfx_res(ki,kj) 
    810         WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ki,kj) 
    811         WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
    812         WRITE(numout,*) ' fheat_res : ', fheat_res(ki,kj) 
    813         WRITE(numout,*) 
    814         WRITE(numout,*) ' - Momentum fluxes ' 
    815         WRITE(numout,*) ' utau      : ', utau(ki,kj)  
    816         WRITE(numout,*) ' vtau      : ', vtau(ki,kj) 
    817       ENDIF 
    818       WRITE(numout,*) ' ' 
    819       ! 
     682      DO ji = mi0(ki), mi1(ki) 
     683         DO jj = mj0(kj), mj1(kj) 
     684 
     685            WRITE(numout,*) ' time step ',kt,' ',cd1             ! print title 
     686 
     687            !---------------- 
     688            !  Simple state 
     689            !---------------- 
     690             
     691            IF ( kn == 1 .OR. kn == -1 ) THEN 
     692               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
     693               WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
     694               WRITE(numout,*) ' Simple state ' 
     695               WRITE(numout,*) ' masks s,u,v   : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 
     696               WRITE(numout,*) ' lat - long    : ', gphit(ji,jj), glamt(ji,jj) 
     697               WRITE(numout,*) ' Time step     : ', numit 
     698               WRITE(numout,*) ' - Ice drift   ' 
     699               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
     700               WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj) 
     701               WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj) 
     702               WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1) 
     703               WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
     704               WRITE(numout,*) ' strength      : ', strength(ji,jj) 
     705               WRITE(numout,*) 
     706               WRITE(numout,*) ' - Cell values ' 
     707               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
     708               WRITE(numout,*) ' cell area     : ', area(ji,jj) 
     709               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
     710               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
     711               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
     712               DO jl = 1, jpl 
     713                  WRITE(numout,*) ' - Category (', jl,')' 
     714                  WRITE(numout,*) ' a_i           : ', a_i(ji,jj,jl) 
     715                  WRITE(numout,*) ' ht_i          : ', ht_i(ji,jj,jl) 
     716                  WRITE(numout,*) ' ht_s          : ', ht_s(ji,jj,jl) 
     717                  WRITE(numout,*) ' v_i           : ', v_i(ji,jj,jl) 
     718                  WRITE(numout,*) ' v_s           : ', v_s(ji,jj,jl) 
     719                  WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)/1.0e9 
     720                  WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)/1.0e9 
     721                  WRITE(numout,*) ' t_su          : ', t_su(ji,jj,jl) 
     722                  WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1,jl) 
     723                  WRITE(numout,*) ' t_i           : ', t_i(ji,jj,1:nlay_i,jl) 
     724                  WRITE(numout,*) ' sm_i          : ', sm_i(ji,jj,jl) 
     725                  WRITE(numout,*) ' smv_i         : ', smv_i(ji,jj,jl) 
     726                  WRITE(numout,*) 
     727               END DO 
     728            ENDIF 
     729            IF( kn == -1 ) THEN 
     730               WRITE(numout,*) ' Mechanical Check ************** ' 
     731               WRITE(numout,*) ' Check what means ice divergence ' 
     732               WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj) 
     733               WRITE(numout,*) ' Total lead fraction     ', ato_i(ji,jj) 
     734               WRITE(numout,*) ' Sum of both             ', ato_i(ji,jj) + at_i(ji,jj) 
     735               WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ji,jj) + at_i(ji,jj) - 1.00 
     736            ENDIF 
     737             
     738 
     739            !-------------------- 
     740            !  Exhaustive state 
     741            !-------------------- 
     742             
     743            IF ( kn .EQ. 2 ) THEN 
     744               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
     745               WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
     746               WRITE(numout,*) ' Exhaustive state ' 
     747               WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 
     748               WRITE(numout,*) ' Time step ', numit 
     749               WRITE(numout,*)  
     750               WRITE(numout,*) ' - Cell values ' 
     751               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
     752               WRITE(numout,*) ' cell area     : ', area(ji,jj) 
     753               WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
     754               WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
     755               WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
     756               WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj) 
     757               WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj) 
     758               WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1) 
     759               WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
     760               WRITE(numout,*) ' strength      : ', strength(ji,jj) 
     761               WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj) 
     762               WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ji,jj)  , ' old_v_ice     : ', old_v_ice(ji,jj)   
     763               WRITE(numout,*) 
     764                
     765               DO jl = 1, jpl 
     766                  WRITE(numout,*) ' - Category (',jl,')' 
     767                  WRITE(numout,*) '   ~~~~~~~~         '  
     768                  WRITE(numout,*) ' ht_i       : ', ht_i(ji,jj,jl)             , ' ht_s       : ', ht_s(ji,jj,jl) 
     769                  WRITE(numout,*) ' t_i        : ', t_i(ji,jj,1:nlay_i,jl) 
     770                  WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl) 
     771                  WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl) 
     772                  WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' old_a_i    : ', old_a_i(ji,jj,jl)    
     773                  WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
     774                  WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' old_v_i    : ', old_v_i(ji,jj,jl)    
     775                  WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
     776                  WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' old_v_s    : ', old_v_s(ji,jj,jl)   
     777                  WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl) 
     778                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ji,jj,1,jl)/1.0e9  
     779                  WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ji,jj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ji,jj,1,jl)/1.0e9 
     780                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ji,jj,2,jl)/1.0e9   
     781                  WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ji,jj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ji,jj,2,jl)/1.0e9 
     782                  WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' old_e_snow : ', old_e_s(ji,jj,1,jl)  
     783                  WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ji,jj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ji,jj,1,jl) 
     784                  WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' old_smv_i  : ', old_smv_i(ji,jj,jl)    
     785                  WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl)  
     786                  WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' old_oa_i   : ', old_oa_i(ji,jj,jl) 
     787                  WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl) 
     788               END DO !jl 
     789                
     790               WRITE(numout,*) 
     791               WRITE(numout,*) ' - Heat / FW fluxes ' 
     792               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
     793               WRITE(numout,*) ' emp        : ', emp      (ji,jj) 
     794               WRITE(numout,*) ' sfx        : ', sfx      (ji,jj) 
     795               WRITE(numout,*) ' sfx_thd    : ', sfx_thd(ji,jj) 
     796               WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ji,jj) 
     797               WRITE(numout,*) ' sfx_mec    : ', sfx_mec  (ji,jj) 
     798               WRITE(numout,*) ' sfx_res    : ', sfx_res(ji,jj) 
     799               WRITE(numout,*) ' fmmec      : ', fmmec    (ji,jj) 
     800               WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj) 
     801               WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj) 
     802               WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj) 
     803               WRITE(numout,*)  
     804               WRITE(numout,*) ' sst        : ', sst_m(ji,jj)   
     805               WRITE(numout,*) ' sss        : ', sss_m(ji,jj)   
     806               WRITE(numout,*)  
     807               WRITE(numout,*) ' - Stresses ' 
     808               WRITE(numout,*) '   ~~~~~~~~ ' 
     809               WRITE(numout,*) ' utau_ice   : ', utau_ice(ji,jj)  
     810               WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ji,jj) 
     811               WRITE(numout,*) ' utau       : ', utau    (ji,jj)  
     812               WRITE(numout,*) ' vtau       : ', vtau    (ji,jj) 
     813               WRITE(numout,*) ' oc. vel. u : ', u_oce   (ji,jj) 
     814               WRITE(numout,*) ' oc. vel. v : ', v_oce   (ji,jj) 
     815            ENDIF 
     816             
     817            !--------------------- 
     818            ! Salt / heat fluxes 
     819            !--------------------- 
     820             
     821            IF ( kn .EQ. 3 ) THEN 
     822               WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
     823               WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
     824               WRITE(numout,*) ' - Salt / Heat Fluxes ' 
     825               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
     826               WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 
     827               WRITE(numout,*) ' Time step ', numit 
     828               WRITE(numout,*) 
     829               WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
     830               WRITE(numout,*) ' qsr       : ', qsr(ji,jj) 
     831               WRITE(numout,*) ' qns       : ', qns(ji,jj) 
     832               WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
     833               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) * r1_rdtice 
     834               WRITE(numout,*) ' qldif     : ', qldif(ji,jj) * r1_rdtice 
     835               WRITE(numout,*) 
     836               WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
     837               WRITE(numout,*) ' emp       : ', emp    (ji,jj) 
     838               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
     839               WRITE(numout,*) ' sfx       : ', sfx    (ji,jj) 
     840               WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj) 
     841               WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ji,jj) 
     842               WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
     843               WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
     844               WRITE(numout,*) 
     845               WRITE(numout,*) ' - Momentum fluxes ' 
     846               WRITE(numout,*) ' utau      : ', utau(ji,jj)  
     847               WRITE(numout,*) ' vtau      : ', vtau(ji,jj) 
     848            ENDIF 
     849            WRITE(numout,*) ' ' 
     850            ! 
     851         END DO 
     852      END DO 
     853 
    820854   END SUBROUTINE lim_prt_state 
    821855 

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

    r4220 r4332  
    335335      !                                                           ! (update freshwater fluxes) 
    336336!RBbug do not understand why see ticket 667 
    337       CALL lbc_lnk( emp, 'T', 1. ) 
     337      !clem-bugsal CALL lbc_lnk( emp, 'T', 1. ) 
    338338      ! 
    339339      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r4220 r4332  
    3333   USE timing         ! Timing 
    3434   USE sbc_ice, ONLY : lk_lim3 
    35  
    3635 
    3736   IMPLICIT NONE 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/step.F90

    r4220 r4332  
    8989      ! Update data, open boundaries, surface boundary condition (including sea-ice) 
    9090      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    91                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
    92  
    9391      IF( lk_tide.AND.(kstp /= nit000 ))   CALL tide_init ( kstp ) 
    9492      IF( lk_tide    )   CALL sbc_tide( kstp ) 
    … …  
    9795      IF( lk_bdy     )   CALL bdy_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries 
    9896 
     97                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
     98                                                      ! clem: moved here for bdy ice purpose 
    9999      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    100100      !  Ocean dynamics : ssh, wn, hdiv, rot                                 ! 

  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/SETTE/sette.sh

    r3708 r4332  
    131131#- 
    132132# Compiler among those in NEMOGCM/ARCH 
    133 COMPILER=PW6_VARGAS 
     133COMPILER=x3750_ADA 
    134134export BATCH_COMMAND_PAR="llsubmit" 
    135135export BATCH_COMMAND_SEQ=$BATCH_COMMAND_PAR 
Note: See TracChangeset for help on using the changeset viewer.