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

Changeset 5051


Ignore:
Timestamp:
2015-02-02T18:31:34+01:00 (9 years ago)
Author:
clem
Message:

LIM3 initialization is now called at the same time as other sbc fields

Location:
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO
Files:
1 added
1 deleted
22 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r5048 r5051  
    1818   PRIVATE 
    1919 
    20    PUBLIC    ice_alloc  !  Called in iceini.F90 
     20   PUBLIC    ice_alloc  !  Called in sbc_lim_init 
    2121 
    2222   !!====================================================================== 
    … …  
    303303   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_melt !: net melting (m) 
    304304 
    305    ! temporary arrays for dummy version of the code 
    306    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dh_i_surf2D, dh_i_bott2D, q_s 
    307  
    308305   !!-------------------------------------------------------------------------- 
    309306   !! * Ice global state variables 
    … …  
    362359   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                      !: ice temperatures 
    363360   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b           !: ice velocity 
    364        
    365  
    366    !!-------------------------------------------------------------------------- 
    367    !! * Increment of global variables 
    368    !!-------------------------------------------------------------------------- 
    369    ! thd refers to changes induced by thermodynamics 
    370    ! trp   ''         ''     ''       advection (transport of ice) 
    371    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_a_i_thd  , d_a_i_trp                 !: icefractions                   
    372    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_v_s_thd  , d_v_s_trp                 !: snow volume 
    373    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_v_i_thd  , d_v_i_trp                 !: ice  volume 
    374    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_smv_i_thd, d_smv_i_trp               !:      
    375    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_sm_i_fl  , d_sm_i_gd                 !: 
    376    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_sm_i_se  , d_sm_i_si  , d_sm_i_la    !: 
    377    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_oa_i_thd , d_oa_i_trp                !: 
    378  
    379    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   d_e_s_thd  , d_e_s_trp     !: 
    380    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   d_e_i_thd  , d_e_i_trp     !: 
    381    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   d_u_ice_dyn, d_v_ice_dyn   !: ice velocity  
    382        
     361             
    383362   !!-------------------------------------------------------------------------- 
    384363   !! * Ice thickness distribution variables 
    … …  
    390369   !! * Ice Run 
    391370   !!-------------------------------------------------------------------------- 
    392    !                                                  !!: ** Namelist namicerun read in iceini ** 
     371   !                                                  !!: ** Namelist namicerun read in sbc_lim_init ** 
    393372   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input) 
    394373   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output) 
    … …  
    400379   !! * Ice diagnostics 
    401380   !!-------------------------------------------------------------------------- 
    402    !! Check if everything down here is necessary 
    403    LOGICAL , PUBLIC                                      ::   ln_limdiahsb  !: flag for ice diag (T) or not (F) 
    404    LOGICAL , PUBLIC                                      ::   ln_limdiaout  !: flag for ice diag (T) or not (F) 
    405    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dv_dt_thd     !: thermodynamic growth rates  
    406    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_vi   !: transport of ice volume 
    407    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_vs   !: transport of snw volume 
    408    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_ei   !: transport of ice enthalpy (W/m2) 
    409    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_es   !: transport of snw enthalpy (W/m2) 
     381   ! Increment of global variables 
     382   ! thd refers to changes induced by thermodynamics 
     383   ! trp   ''         ''     ''       advection (transport of ice) 
     384   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   d_a_i_thd  , d_a_i_trp                 !: icefractions                   
     385   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   d_v_s_thd  , d_v_s_trp                 !: snow volume 
     386   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   d_v_i_thd  , d_v_i_trp                 !: ice  volume 
     387   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   d_smv_i_thd, d_smv_i_trp               !:      
     388   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   d_oa_i_thd , d_oa_i_trp                !: 
     389 
     390   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   d_e_s_thd  , d_e_s_trp     !: 
     391   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   d_e_i_thd  , d_e_i_trp     !: 
     392   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   d_u_ice_dyn, d_v_ice_dyn   !: ice velocity  
     393 
     394   LOGICAL , PUBLIC                                        ::   ln_limdiahsb  !: flag for ice diag (T) or not (F) 
     395   LOGICAL , PUBLIC                                        ::   ln_limdiaout  !: flag for ice diag (T) or not (F) 
     396   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_vi   !: transport of ice volume 
     397   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_vs   !: transport of snw volume 
     398   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_ei   !: transport of ice enthalpy (W/m2) 
     399   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_es   !: transport of snw enthalpy (W/m2) 
    410400   ! 
    411    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_heat_dhc !: snw/ice heat content variation   [W/m2]  
     401   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_heat_dhc !: snw/ice heat content variation   [W/m2]  
    412402   ! 
    413403   INTEGER , PUBLIC ::   jiindx, jjindx        !: indexes of the debugging point 
    … …  
    502492      ALLOCATE( d_a_i_thd(jpi,jpj,jpl) , d_a_i_trp (jpi,jpj,jpl) , d_v_s_thd  (jpi,jpj,jpl) , d_v_s_trp  (jpi,jpj,jpl) ,   & 
    503493         &      d_v_i_thd(jpi,jpj,jpl) , d_v_i_trp (jpi,jpj,jpl) , d_smv_i_thd(jpi,jpj,jpl) , d_smv_i_trp(jpi,jpj,jpl) ,   &      
    504          &      d_sm_i_fl(jpi,jpj,jpl) , d_sm_i_gd (jpi,jpj,jpl) , d_sm_i_se  (jpi,jpj,jpl) , d_sm_i_si  (jpi,jpj,jpl) ,   & 
    505          &      d_sm_i_la(jpi,jpj,jpl) , d_oa_i_thd(jpi,jpj,jpl) , d_oa_i_trp (jpi,jpj,jpl) ,   & 
     494         &      d_oa_i_thd(jpi,jpj,jpl) , d_oa_i_trp (jpi,jpj,jpl) ,   & 
    506495         &     STAT=ierr(ii) ) 
    507496      ii = ii + 1 
    … …  
    515504      ! * Ice diagnostics 
    516505      ii = ii + 1 
    517       ALLOCATE( dv_dt_thd(jpi,jpj,jpl),    & 
    518          &      diag_trp_vi(jpi,jpj), diag_trp_vs  (jpi,jpj), diag_trp_ei(jpi,jpj),   &  
     506      ALLOCATE( diag_trp_vi(jpi,jpj), diag_trp_vs  (jpi,jpj), diag_trp_ei(jpi,jpj),   &  
    519507         &      diag_trp_es(jpi,jpj), diag_heat_dhc(jpi,jpj),  STAT=ierr(ii) ) 
    520508 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90

    r4333 r5051  
    2828   USE lbclnk 
    2929   USE timing          ! Timing 
     30   USE wrk_nemo 
    3031 
    3132   IMPLICIT NONE 
    … …  
    3738CONTAINS 
    3839 
    39    SUBROUTINE lim_cat_1D(zhti,zhts,zai,zht_i,zht_s,za_i) 
     40   SUBROUTINE lim_cat_1D( zhti, zhts, zai, zht_i, zht_s, za_i ) 
    4041      !! Local variables 
    4142      INTEGER  :: ji, jk, jl             ! dummy loop indices 
    42       INTEGER  :: ijpij, i_fill, jl0, ztest_1, ztest_2, ztest_3, ztest_4, ztests   
    43       REAL(wp) :: zarg, zV, zconv 
     43      INTEGER  :: ijpij, i_fill, jl0   
     44      REAL(wp) :: zarg, zV, zconv, zdh 
    4445      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables 
    4546      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables 
    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 
     47      INTEGER , POINTER, DIMENSION(:)         ::   itest 
    4948  
    5049      IF( nn_timing == 1 )  CALL timing_start('limcat_1D') 
     50 
     51      CALL wrk_alloc( 4, itest ) 
    5152      !-------------------------------------------------------------------- 
    5253      ! initialisation of variables 
    … …  
    6667      !         cannot be filled, so we try to fill jpl-1 categories... 
    6768      !         And so forth iteratively until the number of categories filled 
    68       !         fulfills ice volume concervation between input and output (ztests=4)  
     69      !         fulfills ice volume concervation between input and output (itest=4)  
    6970      !-------------------------------------------------------------------------------------- 
    7071 
    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 
    82   
    8372      ! ---------------------------------------- 
    8473      ! distribution over the jpl ice categories 
    … …  
    8978 
    9079         ! initialisation of tests 
    91          ztest_1 = 0 
    92          ztest_2 = 0 
    93          ztest_3 = 0 
    94          ztest_4 = 0 
    95          ztests  = 0 
     80         itest(:)  = 0 
    9681          
    97          i_fill = jpl + 1                                    !==================================== 
    98          DO WHILE ( ( ztests /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories   
    99             ! iteration                                      !==================================== 
     82         i_fill = jpl + 1                                             !==================================== 
     83         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories   
     84            ! iteration                                               !==================================== 
    10085            i_fill = i_fill - 1 
    10186             
    … …  
    10893               zht_i(ji,1) = zhti(ji) 
    10994               za_i (ji,1) = zai (ji) 
    110                ! *** case ice is thicker: fill categories >1 
     95 
     96            ! *** case ice is thicker: fill categories >1 
    11197            ELSE 
    11298 
    113                ! Fill ice thicknesses except the last one (i_fill) by (hmax-hmin)/2  
     99               ! Fill ice thicknesses except the last one (i_fill) by hmean  
    114100               DO jl = 1, i_fill - 1 
    115                   zht_i(ji,jl) = ( zhi_max(jl) + zhi_max(jl-1) ) * 0.5_wp 
     101                  zht_i(ji,jl) = hi_mean(jl) 
    116102               END DO 
    117103                
    … …  
    119105               jl0 = i_fill 
    120106               DO jl = 1, i_fill 
    121                   IF ( ( zhti(ji) >= zhi_max(jl-1) ) .AND. ( zhti(ji) < zhi_max(jl) ) ) THEN 
     107                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
    122108                     jl0 = jl 
    123109           CYCLE 
    … …  
    129115               DO jl = 1, i_fill - 1 
    130116                  IF ( jl == jl0 ) CYCLE 
    131                   zarg           = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
     117                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
    132118                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
    133119               END DO 
    … …  
    138124               ! Ice thickness in the last (i_fill) category 
    139125               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 
    140                zht_i(ji,i_fill) = ( zhti(ji)*zai(ji) - zV ) / za_i(ji,i_fill)  
     126               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill)  
    141127                
    142128            ENDIF ! case ice is thick or thin 
    … …  
    147133            ! Test 1: area conservation 
    148134            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 
    149             IF ( zconv < epsi06 ) ztest_1 = 1 
     135            IF ( zconv < epsi06 ) itest(1) = 1 
    150136             
    151137            ! Test 2: volume conservation 
    152138            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 
    153             IF ( zconv < epsi06 ) ztest_2 = 1 
     139            IF ( zconv < epsi06 ) itest(2) = 1 
    154140             
    155141            ! Test 3: thickness of the last category is in-bounds ? 
    156             IF ( zht_i(ji,i_fill) >= zhi_max(i_fill-1) ) ztest_3 = 1 
     142            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 
    157143             
    158144            ! Test 4: positivity of ice concentrations 
    159             ztest_4 = 1 
     145            itest(4) = 1 
    160146            DO jl = 1, i_fill 
    161                IF ( za_i(ji,jl) < 0._wp ) ztest_4 = 0 
    162             END DO 
    163              
    164             ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 
     147               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 
     148            END DO             
    165149                                                           !============================ 
    166150         END DO                                            ! end iteration on categories 
    167151                                                           !============================ 
    168          ! Check if tests have passed (i.e. volume conservation...) 
    169          !IF ( ztests /= 4 ) THEN 
    170          !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
    171          !   WRITE(numout,*) ' !! ALERT categories distribution !!' 
    172          !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
    173          !   WRITE(numout,*) ' *** ztests is not equal to 4 ' 
    174          !   WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
    175          !   WRITE(numout,*) 'i_fill=',i_fill 
    176          !   WRITE(numout,*) 'zai(ji)=',zai(ji) 
    177          !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:) 
    178          !ENDIF 
    179           
    180152         ENDIF ! if zhti > 0 
    181  
    182153      END DO ! i loop 
    183154 
    … …  
    193164               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )  
    194165               ! recompute ht_i, ht_s avoiding out of bounds values 
    195                zht_i(ji,jl) = MIN( zhi_max(jl), zht_i(ji,jl) + zdh ) 
     166               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh ) 
    196167               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn ) 
    197168            ENDIF 
    … …  
    199170      ENDDO 
    200171 
     172      CALL wrk_dealloc( 4, itest ) 
     173      ! 
    201174      IF( nn_timing == 1 )  CALL timing_stop('limcat_1D') 
    202175      

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

    r4990 r5051  
    3232 
    3333   PUBLIC   lim_diahsb        ! routine called by ice_step.F90 
    34    !!PUBLIC   lim_diahsb_init   ! routine called by ice_init.F90 
    35    !!PUBLIC   lim_diahsb_rst   ! routine called by ice_init.F90 
    3634 
    3735   real(wp) ::   frc_sal, frc_vol   ! global forcing trends 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r4990 r5051  
    8787      !! * Local variables 
    8888      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices 
    89       REAL(wp)   :: epsi20, ztmelts, zdh 
     89      REAL(wp)   :: ztmelts, zdh 
    9090      INTEGER    :: i_hemis, i_fill, jl0   
    9191      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 
    … …  
    100100      CALL wrk_alloc( jpl,   2, zh_i_ini,  za_i_ini,  zv_i_ini ) 
    101101      CALL wrk_alloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 
    102  
    103       epsi20   = 1.e-20_wp 
    104102 
    105103      IF(lwp) WRITE(numout,*) 
    … …  
    197195               !--- Ice thicknesses in the i_fill - 1 first categories 
    198196               DO jl = 1, i_fill - 1 
    199                   zh_i_ini(jl,i_hemis)    = 0.5 * ( hi_max(jl) + hi_max(jl-1) ) 
     197                  zh_i_ini(jl,i_hemis) = hi_mean(jl) 
    200198               END DO 
    201199                
    202200               !--- jl0: most likely index where cc will be maximum 
    203201               DO jl = 1, jpl 
    204                   IF ( ( zht_i_ini(i_hemis) .GT. hi_max(jl-1) ) .AND. & 
    205                      ( zht_i_ini(i_hemis) .LE. hi_max(jl)   ) ) THEN 
     202                  IF ( ( zht_i_ini(i_hemis) > hi_max(jl-1) ) .AND. & 
     203                     & ( zht_i_ini(i_hemis) <= hi_max(jl)   ) ) THEN 
    206204                     jl0 = jl 
    207205                  ENDIF 
    … …  
    267265 
    268266            ! Test 3: thickness of the last category is in-bounds ? 
    269             IF ( zh_i_ini(i_fill, i_hemis) .GT. hi_max(i_fill-1) ) THEN 
     267            IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN 
    270268               ztest_3 = 1 
    271269            ELSE 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r4990 r5051  
    2727   USE wrk_nemo         ! work arrays 
    2828   USE prtctl           ! Print control 
    29   ! Check budget (Rousset) 
     29 
    3030   USE iom              ! I/O manager 
    3131   USE lib_fortran      ! glob_sum 
    … …  
    4040   PUBLIC   lim_itd_me_icestrength 
    4141   PUBLIC   lim_itd_me_init 
    42    PUBLIC   lim_itd_me_zapsmall 
    43    PUBLIC   lim_itd_me_alloc        ! called by iceini.F90 
     42   PUBLIC   lim_itd_me_alloc        ! called by sbc_lim_init  
    4443 
    4544   !----------------------------------------------------------------------- 
    … …  
    375374      !-----------------------------------------------------------------------------! 
    376375      CALL lim_var_glo2eqv 
    377       CALL lim_itd_me_zapsmall 
     376      CALL lim_var_zapsmall 
    378377 
    379378 
    … …  
    632631 
    633632      !     ! Zero out categories with very small areas 
    634       CALL lim_itd_me_zapsmall 
     633      CALL lim_var_zapsmall 
    635634 
    636635      !------------------------------------------------------------------------------! 
    … …  
    13751374   END SUBROUTINE lim_itd_me_init 
    13761375 
    1377  
    1378    SUBROUTINE lim_itd_me_zapsmall 
    1379       !!------------------------------------------------------------------- 
    1380       !!                   ***  ROUTINE lim_itd_me_zapsmall *** 
    1381       !! 
    1382       !! ** Purpose :   Remove too small sea ice areas and correct salt fluxes 
    1383       !! 
    1384       !! history : 
    1385       !! author: William H. Lipscomb, LANL 
    1386       !! Nov 2003:  Modified by Julie Schramm to conserve volume and energy 
    1387       !! Sept 2004: Modified by William Lipscomb; replaced normalize_state with 
    1388       !!            additions to local freshwater, salt, and heat fluxes 
    1389       !!  9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code 
    1390       !!------------------------------------------------------------------- 
    1391       INTEGER ::   ji, jj, jl, jk   ! dummy loop indices 
    1392       INTEGER ::   icells           ! number of cells with ice to zap 
    1393  
    1394       REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace 
    1395       REAL(wp)                          ::   zmask_glo, zsal, zvi, zvs, zei, zes 
    1396 !!gm      REAL(wp) ::   xtmp      ! temporary variable 
    1397       !!------------------------------------------------------------------- 
    1398  
    1399       CALL wrk_alloc( jpi, jpj, zmask ) 
    1400  
    1401       ! to be sure that at_i is the sum of a_i(jl) 
    1402       at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    1403  
    1404       DO jl = 1, jpl 
    1405          !----------------------------------------------------------------- 
    1406          ! Count categories to be zapped. 
    1407          !----------------------------------------------------------------- 
    1408          icells = 0 
    1409          zmask(:,:)  = 0._wp 
    1410          DO jj = 1, jpj 
    1411             DO ji = 1, jpi 
    1412                IF( a_i(ji,jj,jl) <= epsi10 .OR. v_i(ji,jj,jl) <= epsi10 .OR. at_i(ji,jj) <= epsi10 ) THEN 
    1413                   zmask(ji,jj) = 1._wp 
    1414                ENDIF 
    1415             END DO 
    1416          END DO 
    1417          !zmask_glo = glob_sum(zmask) 
    1418          !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean ' 
    1419  
    1420          !----------------------------------------------------------------- 
    1421          ! Zap ice energy and use ocean heat to melt ice 
    1422          !----------------------------------------------------------------- 
    1423  
    1424          DO jk = 1, nlay_i 
    1425             DO jj = 1 , jpj 
    1426                DO ji = 1 , jpi 
    1427                   zei  = e_i(ji,jj,jk,jl) 
    1428                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) 
    1429                   t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) + rtt * zmask(ji,jj) 
    1430                   ! update exchanges with ocean 
    1431                   hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    1432                END DO 
    1433             END DO 
    1434          END DO 
    1435  
    1436          DO jj = 1 , jpj 
    1437             DO ji = 1 , jpi 
    1438                 
    1439                zsal = smv_i(ji,jj,jl) 
    1440                zvi  = v_i(ji,jj,jl) 
    1441                zvs  = v_s(ji,jj,jl) 
    1442                zes  = e_s(ji,jj,1,jl) 
    1443                !----------------------------------------------------------------- 
    1444                ! Zap snow energy and use ocean heat to melt snow 
    1445                !----------------------------------------------------------------- 
    1446                !           xtmp = esnon(i,j,n) / dt ! < 0 
    1447                !           fhnet(i,j)      = fhnet(i,j)      + xtmp 
    1448                !           fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp 
    1449                ! xtmp is greater than 0 
    1450                ! fluxes are positive to the ocean 
    1451                ! here the flux has to be negative for the ocean 
    1452                t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 
    1453  
    1454                !----------------------------------------------------------------- 
    1455                ! zap ice and snow volume, add water and salt to ocean 
    1456                !----------------------------------------------------------------- 
    1457                ato_i(ji,jj)    = a_i  (ji,jj,jl) *           zmask(ji,jj)   + ato_i(ji,jj) 
    1458                a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1459                v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1460                v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1461                t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 
    1462                oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1463                smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1464                e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 
    1465                ! additional condition 
    1466                IF( v_s(ji,jj,jl) <= epsi10 ) THEN 
    1467                   v_s(ji,jj,jl)   = 0._wp 
    1468                   e_s(ji,jj,1,jl) = 0._wp 
    1469                ENDIF 
    1470                ! update exchanges with ocean 
    1471                sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
    1472                wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
    1473                wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
    1474                hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    1475             END DO 
    1476          END DO 
    1477       END DO ! jl  
    1478  
    1479       ! to be sure that at_i is the sum of a_i(jl) 
    1480       at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    1481       ! 
    1482       CALL wrk_dealloc( jpi, jpj, zmask ) 
    1483       ! 
    1484    END SUBROUTINE lim_itd_me_zapsmall 
    1485  
    14861376#else 
    14871377   !!---------------------------------------------------------------------- 
    … …  
    14971387   SUBROUTINE lim_itd_me_init 
    14981388   END SUBROUTINE lim_itd_me_init 
    1499    SUBROUTINE lim_itd_me_zapsmall 
    1500    END SUBROUTINE lim_itd_me_zapsmall 
    15011389#endif 
    15021390   !!====================================================================== 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r5048 r5051  
    954954         zshiftflag = 0 
    955955 
    956 !clem-change 
    957956         DO jj = 1, jpj 
    958957            DO ji = 1, jpi 
    … …  
    976975            zdvice(:,:,jl) = 0._wp 
    977976         ENDIF 
    978 !clem-change 
    979977 
    980978!         ! clem-change begin: why not doing that? 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90

    r4161 r5051  
    2323   PRIVATE 
    2424 
    25    PUBLIC   lim_msh   ! routine called by ice_ini.F90 
     25   PUBLIC   lim_msh   ! routine called by sbcice_lim.F90 
    2626 
    2727   !!---------------------------------------------------------------------- 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r5047 r5051  
    129129 
    130130      REAL(wp) ::   zresm         ! Maximal error on ice velocity 
    131       REAL(wp) ::   zdummy        ! dummy argument 
    132131      REAL(wp) ::   zintb, zintn  ! dummy argument 
    133132 
    … …  
    636635!CDIR NOVERRCHK 
    637636         DO ji = fs_2, fs_jpim1 
    638             zdummy = vt_i(ji,jj) 
    639             IF ( zdummy .LE. hminrhg ) THEN 
     637            IF ( vt_i(ji,jj) .LE. hminrhg ) THEN 
    640638               u_ice(ji,jj) = u_oce(ji,jj) 
    641639               v_ice(ji,jj) = v_oce(ji,jj) 
    642             ENDIF ! zdummy 
     640            ENDIF 
    643641         END DO 
    644642      END DO 
    … …  
    657655      DO jj = k_j1+1, k_jpj-1  
    658656         DO ji = fs_2, fs_jpim1 
    659             zdummy = vt_i(ji,jj) 
    660             IF ( zdummy .LE. hminrhg ) THEN 
     657            IF ( vt_i(ji,jj) .LE. hminrhg ) THEN 
    661658               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   & 
    662659                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 
    … …  
    666663                  &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 
    667664                  &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
    668             ENDIF ! zdummy 
     665            ENDIF  
    669666         END DO 
    670667      END DO 
    … …  
    680677            !- divu_i(:,:), zdt(:,:): divergence and tension at centre  
    681678            !- zds(:,:): shear on northeast corner of grid cells 
    682             zdummy = vt_i(ji,jj) 
    683             IF ( zdummy .LE. hminrhg ) THEN 
     679            IF ( vt_i(ji,jj) .LE. hminrhg ) THEN 
    684680 
    685681               divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
    … …  
    718714               delta_i(ji,jj) = delta + creepl 
    719715             
    720             ENDIF ! zdummy 
     716            ENDIF 
    721717 
    722718         END DO !jj 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90

    r4990 r5051  
    3333   PUBLIC   lim_rst_opn    ! routine called by icestep.F90 
    3434   PUBLIC   lim_rst_write  ! routine called by icestep.F90 
    35    PUBLIC   lim_rst_read   ! routine called by iceini.F90 
     35   PUBLIC   lim_rst_read   ! routine called by sbc_lim_init 
    3636 
    3737   LOGICAL, PUBLIC ::   lrst_ice         !: logical to control the ice restart write  
    … …  
    165165      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i'    , stress2_i  ) 
    166166      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i'   , stress12_i ) 
    167       CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass )   !clem modif 
    168       CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) !clem modif 
     167      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass ) 
     168      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) 
    169169 
    170170      DO jl = 1, jpl  
    … …  
    395395      CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  ) 
    396396      CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i ) 
    397       CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass )   !clem modif 
    398       CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) !clem modif 
     397      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass ) 
     398      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) 
    399399 
    400400      DO jl = 1, jpl  

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r5020 r5051  
    4040   USE prtctl           ! Print control 
    4141   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    42    USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget 
     42   USE traqsr           ! add penetration of solar flux in the calculation of heat budget 
    4343   USE iom 
    4444   USE domvvl           ! Variable volume 
    … …  
    4747   PRIVATE 
    4848 
    49    PUBLIC   lim_sbc_init   ! called by ice_init 
     49   PUBLIC   lim_sbc_init   ! called by sbc_lim_init 
    5050   PUBLIC   lim_sbc_flx    ! called by sbc_ice_lim 
    5151   PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim 
    … …  
    343343         END WHERE 
    344344      ENDIF 
    345       ! clem modif 
     345       
    346346      IF( .NOT. ln_rstart ) THEN 
    347347         fraqsr_1lev(:,:) = 1._wp 
    348348      ENDIF 
    349349      ! 
    350       ! clem: snwice_mass in the restart file now 
    351350      IF( .NOT. ln_rstart ) THEN 
    352351         !                                      ! embedded sea ice 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5049 r5051  
    4949 
    5050   PUBLIC   lim_thd        ! called by limstp module 
    51    PUBLIC   lim_thd_init   ! called by iceini module 
     51   PUBLIC   lim_thd_init   ! called by sbc_lim_init 
    5252 
    5353   !! * Substitutions 
    … …  
    9292      ! 
    9393      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    94       REAL(wp) :: zda 
    9594      ! 
    9695      REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns 
    … …  
    363362            !-------------------------------- 
    364363 
     364            !--------------------------------------! 
     365            ! --- Ice/Snow Temperature profile --- ! 
     366            !--------------------------------------! 
     367            CALL lim_thd_dif( 1, nbpb ) 
     368 
    365369            !---------------------------------! 
    366             ! Ice/Snow Temperature profile    ! 
    367             !---------------------------------! 
    368             CALL lim_thd_dif( 1, nbpb ) 
    369  
    370             !---------------------------------! 
    371             ! Ice/Snow thicnkess              ! 
     370            ! --- Ice/Snow thickness ---      ! 
    372371            !---------------------------------! 
    373372            CALL lim_thd_dh( 1, nbpb )     
    … …  
    377376                                             
    378377            !---------------------------------! 
    379             ! --- Ice salinity --- ! 
     378            ! --- Ice salinity ---            ! 
    380379            !---------------------------------! 
    381380            CALL lim_thd_sal( 1, nbpb )     
    382381 
    383382            !---------------------------------! 
    384             ! --- temperature update --- ! 
     383            ! --- temperature update ---      ! 
    385384            !---------------------------------! 
    386385            CALL lim_thd_temp( 1, nbpb ) 
     386 
     387            !------------------------------------! 
     388            ! --- lateral melting if monocat --- ! 
     389            !------------------------------------! 
     390            IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN 
     391               CALL lim_thd_lam( 1, nbpb ) 
     392            END IF 
    387393 
    388394            !-------------------------------- 
    … …  
    435441              CALL tab_1d_2d( nbpb, hfx_err       , npb, hfx_err_1d(1:nbpb)   , jpi, jpj ) 
    436442              CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj ) 
    437               CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb)   , jpi, jpj ) 
    438  
    439             IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN 
    440               CALL tab_1d_2d( nbpb, dh_i_melt(:,:,jl) , npb, dh_i_melt_1d(1:nbpb) , jpi, jpj ) 
    441             ENDIF 
     443              CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj ) 
     444 
     445!clem            IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN 
     446!clem              CALL tab_1d_2d( nbpb, dh_i_melt(:,:,jl) , npb, dh_i_melt_1d(1:nbpb) , jpi, jpj ) 
     447!clem            ENDIF 
    442448            ! 
    443449              CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
    … …  
    477483      !---------------------------------- 
    478484      CALL lim_var_eqv2glo 
    479  
    480       !---------------------------------- 
    481       ! 5.X) Lateral melting 
    482       !---------------------------------- 
    483       IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN 
    484  
    485          WRITE(numout,*) ' Lateral melting ON ' 
    486  
    487          ! select points where lateral melting occurs 
    488          jl = 1 
    489  
    490          nbplm = 0 
    491          DO jj = 1, jpj 
    492             DO ji = 1, jpi 
    493                IF ( ( dh_i_melt(ji,jj,jl)                  .LT.-epsi10 ) .AND.    & 
    494      &              ( ht_i(ji,jj,jl) - dh_i_melt(ji,jj,jl) .GT. epsi10 ) .AND.    &  
    495      &              ( ht_i(ji,jj,jl)                       .GT. epsi10 ) ) THEN      
    496                   nbplm   = nbplm  + 1 
    497                   nplm(nbplm) = (jj - 1) * jpi + ji 
    498                ENDIF 
    499             END DO 
    500          END DO 
    501  
    502          IF( nbplm > 0 ) THEN  ! If there is no net melting, do nothing 
    503  
    504             ! Move to 1D arrays 
    505             !------------------------- 
    506  
    507             CALL tab_2d_1d( nbplm, a_i_1d      (1:nbplm), a_i(:,:,jl)       , jpi, jpj, nplm(1:nbplm) ) 
    508             CALL tab_2d_1d( nbplm, ht_i_1d     (1:nbplm), ht_i(:,:,jl)      , jpi, jpj, nplm(1:nbplm) ) 
    509             CALL tab_2d_1d( nbplm, dh_i_melt_1d(1:nbplm), dh_i_melt(:,:,jl) , jpi, jpj, nplm(1:nbplm) ) 
    510  
    511             ! Compute lateral melting (dA = A/2h dh ) 
    512             DO ji = 1, nbplm 
    513                zda        = a_i_1d(ji) * dh_i_melt_1d(ji) / ( 2._wp * ht_i_1d(ji) ) 
    514                a_i_1d(ji) = a_i_1d(ji) + zda 
    515             END DO 
    516  
    517             ! Move back to 2D arrays 
    518             !------------------------- 
    519             CALL tab_1d_2d( nbplm, a_i (:,:,jl)  , nplm, a_i_1d     (1:nbplm)   , jpi, jpj ) 
    520             at_i(:,:) = a_i(:,:,jl) 
    521  
    522          ENDIF 
    523  
    524       ENDIF 
    525485 
    526486      !-------------------------------------------- 
    … …  
    602562   END SUBROUTINE lim_thd_temp 
    603563 
     564   SUBROUTINE lim_thd_lam( kideb, kiut ) 
     565      !!----------------------------------------------------------------------- 
     566      !!                   ***  ROUTINE lim_thd_lam ***  
     567      !!                  
     568      !! ** Purpose :   Lateral melting in case monocategory 
     569      !!                          ( dA = A/2h dh ) 
     570      !!----------------------------------------------------------------------- 
     571      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop 
     572      INTEGER             ::   ji            ! dummy loop indices 
     573 
     574      WRITE(numout,*) ' Lateral melting ON ' 
     575      DO ji = kideb, kiut 
     576         IF( ht_i_1d(ji) > epsi10 .AND. dh_i_melt_1d(ji) < 0._wp ) THEN      
     577            a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * dh_i_melt_1d(ji) / ( 2._wp * ht_i_1d(ji) ) )  
     578         END IF 
     579      END DO 
     580      at_i_1d(:) = a_i_1d(:) 
     581       
     582   END SUBROUTINE lim_thd_lam 
     583 
    604584   SUBROUTINE lim_thd_init 
    605585      !!----------------------------------------------------------------------- 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5048 r5051  
    412412      ENDIF 
    413413 
    414       !clem debug. Just to be sure that enthalpy at nlay_i+1 is null 
     414      ! Just to be sure that enthalpy at nlay_i+1 is null 
    415415      DO ji = kideb, kiut 
    416416         q_i_1d(ji,nlay_i+1) = 0._wp 
    … …  
    694694      !------------------------------------------------------------- 
    695695      ! 
    696       IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN 
     696      IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN 
    697697         DO ji = kideb, kiut 
    698698            dh_i_melt_1d(ji) = MIN( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji), 0._wp ) 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5048 r5051  
    116116      REAL(wp) ::   zhsu 
    117117       
    118       REAL(wp), POINTER, DIMENSION(:)     ::   ztfs        ! ice melting point 
    119118      REAL(wp), POINTER, DIMENSION(:)     ::   ztsub       ! old surface temperature (before the iterative procedure ) 
    120119      REAL(wp), POINTER, DIMENSION(:)     ::   ztsubit     ! surface temperature at previous iteration 
    … …  
    170169      !  
    171170      CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 
    172       CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
     171      CALL wrk_alloc( jpij, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    173172      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 
    174173      CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 
    … …  
    189188      ! 1) Initialization                                                            ! 
    190189      !------------------------------------------------------------------------------! 
    191       ! clem clean: replace just ztfs by rtt 
    192190      DO ji = kideb , kiut 
    193191         ! is there snow or not 
    194192         isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ) 
    195          ! surface temperature of fusion 
    196          ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 
    197193         ! layer thickness 
    198194         zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i ) 
    … …  
    296292         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr. 
    297293         ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter 
    298          t_su_1d   (ji) =  MIN( t_su_1d(ji), ztfs(ji) - ztsu_err )  ! necessary 
    299          zerrit   (ji) =  1000._wp                                ! initial value of error 
     294         t_su_1d   (ji) =  MIN( t_su_1d(ji), rtt - ztsu_err )    ! necessary 
     295         zerrit   (ji) =  1000._wp                               ! initial value of error 
    300296      END DO 
    301297 
    … …  
    723719            isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) )  )  ) 
    724720            ztsubit(ji) = t_su_1d(ji) 
    725             IF( t_su_1d(ji) < ztfs(ji) ) & 
     721            IF( t_su_1d(ji) < rtt ) & 
    726722               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   & 
    727723               &          + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
    … …  
    735731         ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 
    736732         DO ji = kideb , kiut 
    737             t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , ztfs(ji) ) , 190._wp  ) 
     733            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rtt ) , 190._wp  ) 
    738734            zerrit(ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )      
    739735         END DO 
    … …  
    832828      ! 
    833829      CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 
    834       CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
     830      CALL wrk_dealloc( jpij, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    835831      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 
    836832      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i,   & 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r4990 r5051  
    8787 
    8888      !-------------------------------------------------------------------------- 
    89       !  1) Cumulative integral of old enthalpy * thicnkess and layers interfaces 
     89      !  1) Cumulative integral of old enthalpy * thickness and layers interfaces 
    9090      !-------------------------------------------------------------------------- 
    9191      zqh_cum0(:,0:nlay_i+2) = 0._wp  

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90

    r5047 r5051  
    3030 
    3131   PUBLIC   lim_thd_sal        ! called by limthd module 
    32    PUBLIC   lim_thd_sal_init   ! called by iceini module 
     32   PUBLIC   lim_thd_sal_init   ! called by sbc_lim_init 
    3333 
    3434   !!---------------------------------------------------------------------- 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r5048 r5051  
    2222   USE limadv         ! ice advection 
    2323   USE limhdf         ! ice horizontal diffusion 
     24   USE limvar         ! clem for ice thickness correction 
     25   ! 
    2426   USE in_out_manager ! I/O manager 
    2527   USE lbclnk         ! lateral boundary conditions -- MPP exchanges 
    … …  
    2830   USE prtctl         ! Print control 
    2931   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    30    USE limvar          ! clem for ice thickness correction 
    31    USE timing          ! Timing 
     32   USE timing         ! Timing 
    3233   USE limcons        ! conservation tests 
    3334 
    … …  
    3536   PRIVATE 
    3637 
    37    PUBLIC   lim_trp    ! called by ice_step 
     38   PUBLIC   lim_trp    ! called by sbcice_lim 
     39 
     40   INTEGER  ::   ncfl                 ! number of ice time step with CFL>1/2   
    3841 
    3942   !! * Substitution 
    … …  
    6063      INTEGER, INTENT(in) ::   kt   ! number of iteration 
    6164      ! 
    62       INTEGER  ::   ji, jj, jk, jl, jn      ! dummy loop indices 
     65      INTEGER  ::   ji, jj, jk, jl, jt   ! dummy loop indices 
    6366      INTEGER  ::   initad                  ! number of sub-timestep for the advection 
    6467      REAL(wp) ::   zcfl , zusnit           !   -      - 
     68      CHARACTER(len=80) :: cltmp 
    6569      ! 
    66       REAL(wp), POINTER, DIMENSION(:,:)      ::   zui_u, zvi_v, zsm, zs0at, zs0ow 
    67       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 
     70      REAL(wp), POINTER, DIMENSION(:,:)      ::   zsm, zs0at 
     71      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e 
     72      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ow 
    6873      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e 
    69       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold   ! old ice volume... 
    70       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zaiold, zhimax   ! old ice concentration and thickness 
    71       REAL(wp), POINTER, DIMENSION(:,:)      ::   zeiold, zesold   ! old enthalpies 
    72       REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei 
     74      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zviold, zvsold           ! old ice volume... 
     75      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness 
     76      REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold   ! old concentration, enthalpies 
     77      REAL(wp) :: zdv, zvi, zvs, zsmv, zes, zei 
    7378      ! 
    74       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
     79      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
    7580      !!--------------------------------------------------------------------- 
    7681      IF( nn_timing == 1 )  CALL timing_start('limtrp') 
    7782 
    78       CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 
    79       CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 
    80       CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, zs0e ) 
    81  
    82       CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold )   ! clem 
     83      CALL wrk_alloc( jpi,jpj,           zsm, zs0at, zatold, zeiold, zesold ) 
     84      CALL wrk_alloc( jpi,jpj,jpl,       zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e ) 
     85      CALL wrk_alloc( jpi,jpj,1,         zs0ow ) 
     86      CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, zs0e ) 
     87      CALL wrk_alloc( jpi,jpj,jpl,       zhimax, zviold, zvsold ) 
    8388 
    8489      IF( numit == nstart .AND. lwp ) THEN 
    … …  
    8893         ENDIF 
    8994         WRITE(numout,*) '~~~~~~~~~~~~' 
     95         ncfl = 0                ! nb of time step with CFL > 1/2 
    9096      ENDIF 
     97 
     98      zsm(:,:) = area(:,:) 
    9199       
    92       zsm(:,:) = area(:,:) 
    93  
    94100      !                             !-------------------------------------! 
    95101      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   ! 
    … …  
    97103 
    98104         ! conservation test 
    99          IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    100  
    101          ! mass and salt flux init (clem) 
    102          zviold(:,:,:)  = v_i(:,:,:) 
    103          zeiold(:,:)  = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )  
    104          zesold(:,:)  = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )  
    105  
    106          !--- Thickness correction init. (clem) ------------------------------- 
     105         IF( ln_limdiahsb )   CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     106 
     107         ! mass and salt flux init 
     108         zviold(:,:,:) = v_i(:,:,:) 
     109         zeiold(:,:)   = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )  
     110         zesold(:,:)   = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )  
     111 
     112         !--- Thickness correction init. ------------------------------- 
    107113         CALL lim_var_glo2eqv 
    108          zaiold(:,:,:) = a_i(:,:,:) 
     114         zatold(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    109115         !--------------------------------------------------------------------- 
    110116         ! Record max of the surrounding ice thicknesses for correction in limupdate 
    … …  
    125131         END DO 
    126132          
    127          !------------------------- 
    128          ! transported fields                                         
    129          !------------------------- 
    130          ! Snow vol, ice vol, salt and age contents, area 
    131          zs0ow(:,:) = ato_i(:,:) * area(:,:)               ! Open water area  
    132          DO jl = 1, jpl 
    133             zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume 
    134             zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume 
    135             zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area 
    136             zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content 
    137             zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content 
    138             zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content 
    139             zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content 
    140          END DO 
    141  
    142          !-------------------------- 
    143          ! Advection of Ice fields  (Prather scheme)                                             
    144          !-------------------------- 
     133         !=============================! 
     134         !==      Prather scheme     ==! 
     135         !=============================! 
     136 
    145137         ! If ice drift field is too fast, use an appropriate time step for advection.          
    146          ! CFL test for stability 
    147          zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) ) 
     138         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )         ! CFL test for stability 
    148139         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) ) 
    149140         IF(lk_mpp )   CALL mpp_max( zcfl ) 
    … …  
    155146         initad = 1 + NINT( MAX( 0._wp, SIGN( 1._wp, zcfl-0.5 ) ) ) 
    156147         zusnit = 1.0 / REAL( initad )  
    157          IF( zcfl > 0.5 .AND. lwp )   & 
    158             WRITE(numout,*) 'lim_trp   : CFL violation at day ', nday, ', cfl = ', zcfl,   & 
    159                &                        ': the ice time stepping is split in two' 
     148         IF( zcfl > 0.5_wp .AND. lwp )   ncfl = ncfl + 1 
     149         IF( numit == nlast .AND. lwp ) THEN 
     150            IF( ncfl > 0 ) THEN    
     151             WRITE(cltmp,'(i6.1)') ncfl 
     152             CALL ctl_stop('STOP',TRIM(cltmp) ) 
     153               CALL ctl_warn( 'lim_trp: ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 
     154            ELSE 
     155               WRITE(numout,*) 'lim_trp : CFL criteria for ice advection is always smaller than 1/2 ' 
     156            ENDIF 
     157         ENDIF 
     158 
     159         !------------------------- 
     160         ! transported fields                                         
     161         !------------------------- 
     162         zs0ow(:,:,1) = ato_i(:,:) * area(:,:)               ! Open water area  
     163         DO jl = 1, jpl 
     164            zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume 
     165            zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume 
     166            zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area 
     167            zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content 
     168            zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content 
     169            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content 
     170            zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content 
     171         END DO 
     172 
    160173 
    161174         IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==! 
    162             DO jn = 1,initad 
    163                CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area 
    164                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    165                CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:),   & 
    166                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     175            DO jt = 1, initad 
     176               CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &             !--- ice open water area 
     177                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     178               CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   & 
     179                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    167180               DO jl = 1, jpl 
    168                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     181                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    169182                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    170183                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   & 
    171184                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    172                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     185                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    173186                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    174187                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   & 
    175188                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    176                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     189                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    177190                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    178191                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   & 
    179192                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    180                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---      
     193                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---      
    181194                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    182195                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   & 
    183196                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    184                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
     197                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations --- 
    185198                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    186199                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &  
    187200                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    188                   CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
     201                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &    !--- snow heat contents --- 
    189202                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    190203                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   & 
    191204                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    192                   DO jk = 1, nlay_i                                                           !--- ice heat contents --- 
    193                      CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     205                  DO jk = 1, nlay_i                                                                !--- ice heat contents --- 
     206                     CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    194207                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    195208                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    … …  
    201214            END DO 
    202215         ELSE 
    203             DO jn = 1, initad 
    204                CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area 
    205                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    206                CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:), sxopw(:,:),   & 
    207                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     216            DO jt = 1, initad 
     217               CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   &             !--- ice open water area 
     218                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
     219               CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:),   & 
     220                  &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  ) 
    208221               DO jl = 1, jpl 
    209                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
     222                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   &    !--- ice volume  --- 
    210223                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    211224                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl),   & 
    212225                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  ) 
    213                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
     226                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  --- 
    214227                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    215228                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl),   & 
    216229                     &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  ) 
    217                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
     230                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity --- 
    218231                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    219232                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl),   & 
    220233                     &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  ) 
    221234 
    222                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
     235                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   &   !--- ice age      --- 
    223236                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    224237                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl),   & 
    225238                     &                                       sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  ) 
    226                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
     239                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations --- 
    227240                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    228241                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a  (:,:,jl), sxa  (:,:,jl),   & 
    229242                     &                                       sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  ) 
    230                   CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
     243                  CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &  !--- snow heat contents --- 
    231244                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    232245                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   & 
    233246                     &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  ) 
    234247                  DO jk = 1, nlay_i                                                           !--- ice heat contents --- 
    235                      CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
     248                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   &  
    236249                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   & 
    237250                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    … …  
    247260         ! Recover the properties from their contents 
    248261         !------------------------------------------- 
    249          zs0ow(:,:) = zs0ow(:,:) / area(:,:) 
    250          DO jl = 1, jpl 
    251             zs0ice(:,:,jl) = zs0ice(:,:,jl) / area(:,:) 
    252             zs0sn (:,:,jl) = zs0sn (:,:,jl) / area(:,:) 
    253             zs0sm (:,:,jl) = zs0sm (:,:,jl) / area(:,:) 
    254             zs0oi (:,:,jl) = zs0oi (:,:,jl) / area(:,:) 
    255             zs0a  (:,:,jl) = zs0a  (:,:,jl) / area(:,:) 
    256             ! 
     262         ato_i(:,:) = zs0ow(:,:,1) / area(:,:) 
     263         DO jl = 1, jpl 
     264            v_i  (:,:,jl)   = zs0ice(:,:,jl) / area(:,:) 
     265            v_s  (:,:,jl)   = zs0sn (:,:,jl) / area(:,:) 
     266            smv_i(:,:,jl)   = zs0sm (:,:,jl) / area(:,:) 
     267            oa_i (:,:,jl)   = zs0oi (:,:,jl) / area(:,:) 
     268            a_i  (:,:,jl)   = zs0a  (:,:,jl) / area(:,:) 
     269            e_s  (:,:,1,jl) = zs0c0 (:,:,jl)               
     270            e_i  (:,:,:,jl) = zs0e  (:,:,:,jl)            
     271         END DO 
     272 
     273         at_i(:,:) = a_i(:,:,1)      ! total ice fraction 
     274         DO jl = 2, jpl 
     275            at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    257276         END DO 
    258277 
    … …  
    261280         !------------------------------------------------------------------------------! 
    262281 
     282         ! 
    263283         !-------------------------------- 
    264284         !  diffusion of open water area 
    265285         !-------------------------------- 
    266          zs0at(:,:) = zs0a(:,:,1)      ! total ice fraction 
    267          DO jl = 2, jpl 
    268             zs0at(:,:) = zs0at(:,:) + zs0a(:,:,jl) 
    269          END DO 
    270          ! 
    271286         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points 
    272287         DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row 
    273288            DO ji = 1 , fs_jpim1   ! vector opt. 
    274                pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji  ,jj) ) ) )   & 
    275                   &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji+1,jj) ) ) ) * ahiu(ji,jj) 
    276                pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0at(ji,jj  ) ) ) )   & 
    277                   &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0at(ji,jj+1) ) ) ) * ahiv(ji,jj) 
     289               pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   & 
     290                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj) 
     291               pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   & 
     292                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj) 
    278293            END DO 
    279294         END DO 
    280295         ! 
    281          CALL lim_hdf( zs0ow (:,:) )   ! Diffusion 
     296         CALL lim_hdf( ato_i (:,:) )   ! Diffusion 
    282297 
    283298         !------------------------------------ 
    … …  
    288303            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row 
    289304               DO ji = 1 , fs_jpim1   ! vector opt. 
    290                   pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji  ,jj,jl) ) ) )   & 
    291                      &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 
    292                   pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -zs0a(ji,jj  ,jl) ) ) )   & 
    293                      &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- zs0a(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 
    294                END DO 
    295             END DO 
    296  
    297             CALL lim_hdf( zs0ice (:,:,jl) ) 
    298             CALL lim_hdf( zs0sn  (:,:,jl) ) 
    299             CALL lim_hdf( zs0sm  (:,:,jl) ) 
    300             CALL lim_hdf( zs0oi  (:,:,jl) ) 
    301             CALL lim_hdf( zs0a   (:,:,jl) ) 
    302             CALL lim_hdf( zs0c0  (:,:,jl) ) 
     305                  pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   & 
     306                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 
     307                  pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   & 
     308                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 
     309               END DO 
     310            END DO 
     311 
     312            CALL lim_hdf( v_i  (:,:,  jl) ) 
     313            CALL lim_hdf( v_s  (:,:,  jl) ) 
     314            CALL lim_hdf( smv_i(:,:,  jl) ) 
     315            CALL lim_hdf( oa_i (:,:,  jl) ) 
     316            CALL lim_hdf( a_i  (:,:,  jl) ) 
     317            CALL lim_hdf( e_s  (:,:,1,jl) ) 
    303318            DO jk = 1, nlay_i 
    304                CALL lim_hdf( zs0e (:,:,jk,jl) ) 
     319               CALL lim_hdf( e_i(:,:,jk,jl) ) 
    305320            END DO 
    306321         END DO 
    … …  
    310325         !------------------------------------------------------------------------------! 
    311326 
     327!!gm & cr   :  MAX should not be active if adv scheme is positive ! 
    312328         !-------------------------------------------------- 
    313329         ! 5.1) Recover mean values over the grid squares. 
    314330         !-------------------------------------------------- 
    315          zs0at(:,:) = 0._wp 
    316331         DO jl = 1, jpl 
    317332            DO jj = 1, jpj 
    318333               DO ji = 1, jpi 
    319                   zs0sn (ji,jj,jl) = MAX( 0._wp, zs0sn (ji,jj,jl) ) 
    320                   zs0ice(ji,jj,jl) = MAX( 0._wp, zs0ice(ji,jj,jl) ) 
    321                   zs0sm (ji,jj,jl) = MAX( 0._wp, zs0sm (ji,jj,jl) ) 
    322                   zs0oi (ji,jj,jl) = MAX( 0._wp, zs0oi (ji,jj,jl) ) 
    323                   zs0a  (ji,jj,jl) = MAX( 0._wp, zs0a  (ji,jj,jl) ) 
    324                   zs0c0 (ji,jj,jl) = MAX( 0._wp, zs0c0 (ji,jj,jl) ) 
    325                   zs0at (ji,jj)    = zs0at(ji,jj) + zs0a(ji,jj,jl) 
    326                END DO 
    327             END DO 
    328          END DO 
    329  
    330          !--------------------------------------------------------- 
    331          ! 5.2) Update and mask variables 
    332          !--------------------------------------------------------- 
    333          DO jl = 1, jpl           
    334             DO jj = 1, jpj 
    335                DO ji = 1, jpi 
    336                   rswitch = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 
    337  
    338                   zvi  = zs0ice(ji,jj,jl) 
    339                   zvs  = zs0sn (ji,jj,jl) 
    340                   zes  = zs0c0 (ji,jj,jl)       
    341                   zsmv = zs0sm (ji,jj,jl) 
    342                   ! 
    343                   ! Remove very small areas 
    344                   v_s(ji,jj,jl)   = rswitch * zs0sn (ji,jj,jl)  
    345                   v_i(ji,jj,jl)   = rswitch * zs0ice(ji,jj,jl) 
    346                   a_i(ji,jj,jl)   = rswitch * zs0a  (ji,jj,jl) 
    347                   e_s(ji,jj,1,jl) = rswitch * zs0c0 (ji,jj,jl)       
    348                   ! Ice salinity and age 
    349                   IF(  num_sal == 2  ) THEN 
    350                      smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), zsmv ), s_i_min * v_i(ji,jj,jl) ) 
    351                   ENDIF 
    352                   oa_i(ji,jj,jl) = MAX( rswitch * zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ), 0._wp ) * a_i(ji,jj,jl) 
    353  
    354                  ! Update fluxes 
    355                   wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice  
    356                   wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
    357                   sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice  
    358                   hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    359               END DO 
    360             END DO 
    361          END DO 
    362  
     334                  v_s  (ji,jj,jl)   = MAX( 0._wp, v_s  (ji,jj,jl) ) 
     335                  v_i  (ji,jj,jl)   = MAX( 0._wp, v_i  (ji,jj,jl) ) 
     336                  smv_i(ji,jj,jl)   = MAX( 0._wp, smv_i(ji,jj,jl) ) 
     337                  oa_i (ji,jj,jl)   = MAX( 0._wp, oa_i (ji,jj,jl) ) 
     338                  a_i  (ji,jj,jl)   = MAX( 0._wp, a_i  (ji,jj,jl) ) 
     339                  e_s  (ji,jj,1,jl) = MAX( 0._wp, e_s  (ji,jj,1,jl) ) 
     340               END DO 
     341            END DO 
     342         END DO 
    363343         DO jl = 1, jpl 
    364344            DO jk = 1, nlay_i 
    365345               DO jj = 1, jpj 
    366346                  DO ji = 1, jpi 
    367                      rswitch          = MAX( 0._wp , SIGN( 1._wp, zs0a(ji,jj,jl) - epsi10 ) ) 
    368                      zei              = zs0e(ji,jj,jk,jl)       
    369                      e_i(ji,jj,jk,jl) = rswitch * MAX( 0._wp, zs0e(ji,jj,jk,jl) ) 
    370                      ! Update fluxes 
    371                      hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    372                   END DO !ji 
    373                END DO ! jj 
    374             END DO ! jk 
    375          END DO ! jl 
    376  
    377          !--- Thickness correction in case too high (clem) -------------------------------------------------------- 
     347                     e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) ) 
     348                  END DO 
     349               END DO 
     350            END DO 
     351         END DO 
     352!!gm & cr  
     353 
     354         ! zap small areas 
     355         CALL lim_var_zapsmall 
     356 
     357         !--- Thickness correction in case too high -------------------------------------------------------- 
    378358         CALL lim_var_glo2eqv 
    379359         DO jl = 1, jpl 
    … …  
    388368                     zei  = SUM( e_i(ji,jj,1:nlay_i,jl) ) 
    389369                     zdv  = v_i(ji,jj,jl) - zviold(ji,jj,jl)    
    390                      !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl)    
    391370                      
    392371                     rswitch = 1._wp 
    393                      IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 
     372                     IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 
    394373                        & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                           
    395374                        ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 
    … …  
    416395                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    417396                  ENDIF 
     397 
    418398               END DO 
    419399            END DO 
    … …  
    466446               ! open water = 1 if at_i=0 
    467447               rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 
    468                ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * zs0ow(ji,jj) 
     448               ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 
    469449            END DO 
    470450         END DO       
    … …  
    507487      ENDIF 
    508488      ! 
    509       CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 
    510       CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 
    511       CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, zs0e ) 
    512  
    513       CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax )   ! clem 
     489      CALL wrk_dealloc( jpi,jpj,           zsm, zs0at, zatold, zeiold, zesold ) 
     490      CALL wrk_dealloc( jpi,jpj,jpl,       zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e ) 
     491      CALL wrk_dealloc( jpi,jpj,1,         zs0ow ) 
     492      CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, zs0e ) 
     493      CALL wrk_dealloc( jpi,jpj,jpl,       zviold, zvsold, zhimax ) 
    514494      ! 
    515495      IF( nn_timing == 1 )  CALL timing_stop('limtrp') 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r4990 r5051  
    8383      ! zap small values 
    8484      !----------------- 
    85       CALL lim_itd_me_zapsmall 
     85      CALL lim_var_zapsmall 
    8686 
    8787      CALL lim_var_glo2eqv 
    … …  
    124124      ! zap small values 
    125125      !----------------- 
    126       CALL lim_itd_me_zapsmall 
     126      CALL lim_var_zapsmall 
    127127 
    128128      !--------------------- 
    … …  
    148148      ! Diagnostics 
    149149      ! ------------------------------------------------- 
     150      DO jl  = 1, jpl 
     151         afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 
     152      END DO 
     153 
    150154      d_u_ice_dyn(:,:)     = u_ice(:,:)     - u_ice_b(:,:) 
    151155      d_v_ice_dyn(:,:)     = v_ice(:,:)     - v_ice_b(:,:) 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r4990 r5051  
    7878      ! zap small values 
    7979      !----------------- 
    80       CALL lim_itd_me_zapsmall 
    81  
     80      CALL lim_var_agg( 1 ) 
     81      CALL lim_var_zapsmall 
    8282      CALL lim_var_glo2eqv 
    8383 
    … …  
    133133      ! zap small values 
    134134      !----------------- 
    135       CALL lim_itd_me_zapsmall 
     135      CALL lim_var_zapsmall 
    136136 
    137137      !--------------------- 
    … …  
    176176      v_ice(:,:) = v_ice(:,:) * tmv(:,:) 
    177177  
     178      ! for outputs 
     179      CALL lim_var_glo2eqv            ! equivalent variables (outputs) 
     180      CALL lim_var_agg(2)             ! aggregate ice thickness categories 
     181 
    178182      ! ------------------------------------------------- 
    179183      ! Diagnostics 
    180184      ! ------------------------------------------------- 
     185      DO jl  = 1, jpl 
     186         afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 
     187      END DO 
     188      afx_tot = afx_thd + afx_dyn 
     189 
    181190      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - a_i_b(:,:,:)  
    182191      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - v_s_b(:,:,:) 
    … …  
    187196      d_smv_i_thd(:,:,:) = 0._wp 
    188197      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:) 
    189       ! diag only (clem) 
    190       dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    191198 
    192199      ! heat content variation (W.m-2) 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r5047 r5051  
    6565   PUBLIC   lim_var_bv           ! 
    6666   PUBLIC   lim_var_salprof1d    ! 
     67   PUBLIC   lim_var_zapsmall 
    6768 
    6869   !!---------------------------------------------------------------------- 
    … …  
    529530   END SUBROUTINE lim_var_salprof1d 
    530531 
     532   SUBROUTINE lim_var_zapsmall 
     533      !!------------------------------------------------------------------- 
     534      !!                   ***  ROUTINE lim_var_zapsmall *** 
     535      !! 
     536      !! ** Purpose :   Remove too small sea ice areas and correct salt fluxes 
     537      !! 
     538      !! history : LIM3.5 - 01-2014 (C. Rousset) original code 
     539      !!------------------------------------------------------------------- 
     540      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
     541 
     542      REAL(wp) ::   zsal, zvi, zvs, zei, zes 
     543      !!------------------------------------------------------------------- 
     544 
     545      DO jl = 1, jpl 
     546 
     547         !----------------------------------------------------------------- 
     548         ! Zap ice energy and use ocean heat to melt ice 
     549         !----------------------------------------------------------------- 
     550         DO jk = 1, nlay_i 
     551            DO jj = 1 , jpj 
     552               DO ji = 1 , jpi 
     553                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
     554                  zei              = e_i(ji,jj,jk,jl) 
     555                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch 
     556                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rtt * ( 1._wp - rswitch ) 
     557                  ! update exchanges with ocean 
     558                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
     559               END DO 
     560            END DO 
     561         END DO 
     562 
     563         DO jj = 1 , jpj 
     564            DO ji = 1 , jpi 
     565               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
     566                
     567               zsal = smv_i(ji,jj,  jl) 
     568               zvi  = v_i  (ji,jj,  jl) 
     569               zvs  = v_s  (ji,jj,  jl) 
     570               zes  = e_s  (ji,jj,1,jl) 
     571               !----------------------------------------------------------------- 
     572               ! Zap snow energy  
     573               !----------------------------------------------------------------- 
     574               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rtt * ( 1._wp - rswitch ) 
     575               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch 
     576 
     577               !----------------------------------------------------------------- 
     578               ! zap ice and snow volume, add water and salt to ocean 
     579               !----------------------------------------------------------------- 
     580               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj) 
     581               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch 
     582               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch 
     583               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch 
     584               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch ) 
     585               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch 
     586               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch 
     587 
     588               ! ice salinity must stay in bounds 
     589               IF(  num_sal == 2  ) THEN 
     590                  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) ) 
     591               ENDIF 
     592               ! update exchanges with ocean 
     593               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
     594               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
     595               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
     596               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
     597            END DO 
     598         END DO 
     599      END DO ! jl  
     600 
     601      ! to be sure that at_i is the sum of a_i(jl) 
     602      at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
     603      ! 
     604   END SUBROUTINE lim_var_zapsmall 
     605 
    531606#else 
    532607   !!---------------------------------------------------------------------- 
    … …  
    546621   SUBROUTINE lim_var_salprof1d    ! Emtpy routines 
    547622   END SUBROUTINE lim_var_salprof1d 
     623   SUBROUTINE lim_var_zapsmall 
     624   END SUBROUTINE lim_var_zapsmall 
    548625#endif 
    549626 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r5048 r5051  
    212212      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow) 
    213213       
    214       CALL iom_put( "afxtot"     , afx_tot              )        ! concentration tendency (total) 
    215       CALL iom_put( "afxdyn"     , afx_dyn              )        ! concentration tendency (dynamics) 
    216       CALL iom_put( "afxthd"     , afx_thd              )        ! concentration tendency (thermo) 
    217  
    218       CALL iom_put ('hfxthd', hfx_thd(:,:) )   !   
    219       CALL iom_put ('hfxdyn', hfx_dyn(:,:) )   !   
    220       CALL iom_put ('hfxres', hfx_res(:,:) )   !   
    221       CALL iom_put ('hfxout', hfx_out(:,:) )   !   
    222       CALL iom_put ('hfxin' , hfx_in(:,:) )    !   
    223       CALL iom_put ('hfxsnw', hfx_snw(:,:) )   !   
    224       CALL iom_put ('hfxsub', hfx_sub(:,:) )   !   
    225       CALL iom_put ('hfxerr', hfx_err(:,:) )   !   
    226       CALL iom_put ('hfxerr_rem', hfx_err_rem(:,:) )   !   
    227        
    228       CALL iom_put ('hfxsum', hfx_sum(:,:) )   !   
    229       CALL iom_put ('hfxbom', hfx_bom(:,:) )   !   
    230       CALL iom_put ('hfxbog', hfx_bog(:,:) )   !   
    231       CALL iom_put ('hfxdif', hfx_dif(:,:) )   !   
    232       CALL iom_put ('hfxopw', hfx_opw(:,:) )   !   
    233       CALL iom_put ('hfxtur', fhtur(:,:) * at_i(:,:) )  ! turbulent heat flux at ice base  
    234       CALL iom_put ('hfxdhc', diag_heat_dhc(:,:) )          ! Heat content variation in snow and ice  
    235       CALL iom_put ('hfxspr', hfx_spr(:,:) )          ! Heat content of snow precip  
     214      CALL iom_put( "afxtot"     , afx_tot * rday       )        ! concentration tendency (total) 
     215      CALL iom_put( "afxdyn"     , afx_dyn * rday       )        ! concentration tendency (dynamics) 
     216      CALL iom_put( "afxthd"     , afx_thd * rday       )        ! concentration tendency (thermo) 
     217 
     218      CALL iom_put ('hfxthd'     , hfx_thd(:,:)        )   !   
     219      CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)        )   !   
     220      CALL iom_put ('hfxres'     , hfx_res(:,:)        )   !   
     221      CALL iom_put ('hfxout'     , hfx_out(:,:)        )   !   
     222      CALL iom_put ('hfxin'      , hfx_in(:,:)          )    !   
     223      CALL iom_put ('hfxsnw'     , hfx_snw(:,:)        )   !   
     224      CALL iom_put ('hfxsub'     , hfx_sub(:,:)        )   !   
     225      CALL iom_put ('hfxerr'     , hfx_err(:,:)        )   !   
     226      CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)    )   !   
     227       
     228      CALL iom_put ('hfxsum'     , hfx_sum(:,:)        )   !   
     229      CALL iom_put ('hfxbom'     , hfx_bom(:,:)        )   !   
     230      CALL iom_put ('hfxbog'     , hfx_bog(:,:)        )   !   
     231      CALL iom_put ('hfxdif'     , hfx_dif(:,:)        )   !   
     232      CALL iom_put ('hfxopw'     , hfx_opw(:,:)        )   !   
     233      CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i(:,:) ) ! turbulent heat flux at ice base  
     234      CALL iom_put ('hfxdhc'     , diag_heat_dhc(:,:)   )   ! Heat content variation in snow and ice  
     235      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  
    236236       
    237237      !-------------------------------- 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r5048 r5051  
    1919   !!---------------------------------------------------------------------- 
    2020   !!   sbc_ice_lim  : sea-ice model time-stepping and update ocean sbc over ice-covered area 
    21    !!   lim_ctl       : alerts in case of ice model crash 
    22    !!   lim_prt_state : ice control print at a given grid point 
    2321   !!---------------------------------------------------------------------- 
    2422   USE oce             ! ocean dynamics and tracers 
    … …  
    2624   USE par_ice         ! sea-ice parameters 
    2725   USE ice             ! LIM-3: ice variables 
    28    USE iceini          ! LIM-3: ice initialisation 
     26   USE thd_ice         ! LIM-3: thermodynamical variables 
    2927   USE dom_ice         ! LIM-3: ice domain 
    3028 
    … …  
    4745   USE limwri          ! Ice outputs 
    4846   USE limrst          ! Ice restarts 
    49    USE limupdate1       ! update of global variables 
    50    USE limupdate2       ! update of global variables 
     47   USE limupdate1      ! update of global variables 
     48   USE limupdate2      ! update of global variables 
    5149   USE limvar          ! Ice variables switch 
     50 
     51   USE limmsh          ! LIM mesh 
     52   USE limistate       ! LIM initial state 
     53   USE limthd_sal      ! LIM ice thermodynamics: salinity 
    5254 
    5355   USE c1d             ! 1D vertical configuration 
    … …  
    6062   USE prtctl          ! Print control 
    6163   USE lib_fortran     !  
     64   USE limctl 
    6265 
    6366#if defined key_bdy  
    … …  
    6972 
    7073   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90 
    71    PUBLIC lim_prt_state 
     74   PUBLIC sbc_lim_init ! routine called by sbcmod.F90 
    7275    
    7376   !! * Substitutions 
    … …  
    114117      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
    115118 
    116       IF( kt == nit000 ) THEN 
    117          IF(lwp) WRITE(numout,*) 
    118          IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition'  
    119          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping' 
    120          ! 
    121          CALL ice_init 
    122          ! 
    123          IF( ln_nicep ) THEN      ! control print at a given point 
    124             jiindx = 15    ;   jjindx =  44 
    125             IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    126          ENDIF 
    127       ENDIF 
    128  
    129119      !                                        !----------------------! 
    130120      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  ! 
    … …  
    133123         !                                           !----------------! 
    134124         ! 
    135          u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)                     ! mean surface ocean current at ice velocity point 
    136          v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)                    ! (C-grid dynamics :  U- & V-points as the ocean) 
    137          ! 
    138          t_bo(:,:) = ( eos_fzp( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) )  ! masked sea surface freezing temperature [Kelvin] 
    139          !                                                                                  ! (set to rt0 over land) 
    140          !                                           ! Ice albedo 
    141          CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )       
    142  
     125         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)      ! mean surface ocean current at ice velocity point 
     126         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)      ! (C-grid dynamics :  U- & V-points as the ocean) 
     127          
     128         ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 
     129         t_bo(:,:) = ( eos_fzp( sss_m ) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )   
     130         !                                                                                       
     131         ! Ice albedo 
     132         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
    143133         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
    144134 
    … …  
    153143         END SELECT 
    154144          
    155          !                                           ! Mask sea ice surface temperature 
     145         ! Mask sea ice surface temperature (set to rt0 over land) 
    156146         DO jl = 1, jpl 
    157             t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) ) 
     147            t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 
    158148         END DO 
    159149      
    … …  
    210200         v_ice_b(:,:)     = v_ice(:,:) 
    211201 
    212          ! salt, heat and mass fluxes 
    213          sfx    (:,:) = 0._wp   ; 
    214          sfx_bri(:,:) = 0._wp   ;  
    215          sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
    216          sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
    217          sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
    218          sfx_res(:,:) = 0._wp 
    219  
    220          wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
    221          wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
    222          wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
    223          wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
    224          wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
    225          wfx_spr(:,:) = 0._wp   ;    
    226  
    227          afx_tot(:,:) = at_i(:,:) ;  afx_dyn(:,:) = 0._wp 
    228          afx_thd(:,:) = 0._wp 
    229  
    230          hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
    231          hfx_thd(:,:) = 0._wp   ;    
    232          hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
    233          hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
    234          hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
    235          hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
    236          hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
    237          hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
     202                          CALL sbc_lim_flx0          ! set diag of mass, heat and salt fluxes to 0 
    238203 
    239204                          CALL lim_rst_opn( kt )     ! Open Ice restart file 
    240205         ! 
    241          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
     206         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
    242207         ! ---------------------------------------------- 
    243208         ! ice dynamics and transport (except in 1D case) 
    … …  
    245210         IF( .NOT. lk_c1d ) THEN 
    246211 
    247          IF ( ln_limdyn ) afx_dyn(:,:)     = at_i(:,:) 
    248  
    249212                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics ) 
    250213                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion ) 
    251214                          CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
    252          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
    253          IF( nn_monocat /= 2 ) CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting) 
     215         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
     216         IF( nn_monocat /= 2 )   & 
     217            &             CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting) 
    254218                          CALL lim_var_agg( 1 )  
    255219#if defined key_bdy 
    … …  
    257221                          CALL lim_var_glo2eqv            ! equivalent variables 
    258222                          CALL bdy_ice_lim( kt ) 
    259                           CALL lim_itd_me_zapsmall 
     223                          CALL lim_var_zapsmall 
    260224                          CALL lim_var_agg(1) 
    261          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print 
     225         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print 
    262226#endif 
    263227                          CALL lim_update1 
    … …  
    273237                          oa_i_b (:,:,:)   = oa_i (:,:,:) 
    274238                          smv_i_b(:,:,:)   = smv_i(:,:,:) 
    275  
    276          IF ( ln_limdyn ) afx_dyn(:,:)     = ( at_i(:,:) - afx_dyn(:,:) ) * r1_rdtice 
    277                           afx_thd(:,:)     = at_i(:,:) 
    278239  
    279240         ! ---------------------------------------------- 
    … …  
    282243                          CALL lim_var_glo2eqv            ! equivalent variables 
    283244                          CALL lim_var_agg(1)             ! aggregate ice categories 
     245 
    284246                          ! previous lead fraction and ice volume for flux calculations 
    285247                          pfrld(:,:)   = 1._wp - at_i(:,:) 
    … …  
    301263                          zcoef = rdt_ice /rday           !  Ice natural aging 
    302264                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 
    303          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
    304                           CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  ! 
    305                           CALL lim_var_agg( 1 )           ! requested by limupdate 
     265         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
     266                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion 
    306267                          CALL lim_update2                ! Global variables update 
    307268 
    308                           CALL lim_var_glo2eqv            ! equivalent variables (outputs) 
    309                           CALL lim_var_agg(2)             ! aggregate ice thickness categories 
    310          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
     269         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
    311270         ! 
    312271                          CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes 
    313272         ! 
    314          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
    315          ! 
    316          !                                           ! Diagnostics and outputs  
    317          IF (ln_limdiaout) CALL lim_diahsb 
    318  
    319                           afx_thd(:,:)     = ( at_i(:,:) - afx_thd(:,:) ) * r1_rdtice 
    320                           afx_tot(:,:)     = ( at_i(:,:) - afx_tot(:,:) ) * r1_rdtice 
     273         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
     274         ! 
     275         IF(ln_limdiaout) CALL lim_diahsb                 ! Diagnostics and outputs  
    321276 
    322277                          CALL lim_wri( 1  )              ! Ice outputs  
    … …  
    342297      IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents 
    343298!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    344  
    345299      ! 
    346300      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim') 
    … …  
    348302   END SUBROUTINE sbc_ice_lim 
    349303    
     304 
     305   SUBROUTINE sbc_lim_init 
     306      !!---------------------------------------------------------------------- 
     307      !!                  ***  ROUTINE sbc_lim_init  *** 
     308      !! 
     309      !! ** purpose :   Allocate all the dynamic arrays of the LIM-3 modules 
     310      !!---------------------------------------------------------------------- 
     311      INTEGER :: ierr 
     312      !!---------------------------------------------------------------------- 
     313      IF(lwp) WRITE(numout,*) 
     314      IF(lwp) WRITE(numout,*) 'sbc_ice_lim : update ocean surface boudary condition'  
     315      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping' 
     316      ! 
     317      ! 
     318      !                                ! Allocate the ice arrays 
     319      ierr =        ice_alloc        ()      ! ice variables 
     320      ierr = ierr + dom_ice_alloc    ()      ! domain 
     321      ierr = ierr + sbc_ice_alloc    ()      ! surface forcing 
     322      ierr = ierr + thd_ice_alloc    ()      ! thermodynamics 
     323      ierr = ierr + lim_itd_me_alloc ()      ! ice thickness distribution - mechanics 
     324      ! 
     325      IF( lk_mpp    )   CALL mpp_sum( ierr ) 
     326      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'sbc_lim_init : unable to allocate ice arrays') 
     327      ! 
     328      !                                ! adequation jpk versus ice/snow layers/categories 
     329      IF( jpl > jpk .OR. (nlay_i+1) > jpk .OR. nlay_s > jpk )   & 
     330         &      CALL ctl_stop( 'STOP',                     & 
     331         &     'sbc_lim_init: the 3rd dimension of workspace arrays is too small.',   & 
     332         &     'use more ocean levels or less ice/snow layers/categories.' ) 
     333 
     334                                       ! Open the reference and configuration namelist files and namelist output file  
     335      CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )  
     336      CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg',    'OLD',     'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
     337      IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 ) 
     338      ! 
     339      CALL ice_run                     ! set some ice run parameters 
     340      ! 
     341      CALL lim_thd_init                ! set ice thermodynics parameters 
     342      ! 
     343      CALL lim_thd_sal_init            ! set ice salinity parameters 
     344      ! 
     345      rdt_ice   = nn_fsbc * rdttra(1)  ! sea-ice timestep 
     346      r1_rdtice = 1._wp / rdt_ice      ! sea-ice timestep inverse 
     347      ! 
     348      CALL lim_msh                     ! ice mesh initialization 
     349      ! 
     350      CALL lim_itd_init                 ! ice thickness distribution initialization 
     351      ! 
     352      CALL lim_itd_me_init             ! ice thickness distribution initialization 
     353      !                                ! Initial sea-ice state 
     354      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst 
     355         numit = 0 
     356         numit = nit000 - 1 
     357         CALL lim_istate 
     358         CALL lim_var_agg(1) 
     359         CALL lim_var_glo2eqv 
     360      ELSE                                    ! start from a restart file 
     361         CALL lim_rst_read 
     362         numit = nit000 - 1 
     363         CALL lim_var_agg(1) 
     364         CALL lim_var_glo2eqv 
     365      ENDIF 
     366      ! 
     367      CALL lim_sbc_init                 ! ice surface boundary condition    
     368      ! 
     369      fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction 
     370      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu 
     371      ! 
     372      nstart = numit  + nn_fsbc       
     373      nitrun = nitend - nit000 + 1  
     374      nlast  = numit  + nitrun  
     375      ! 
     376      IF( nstock == 0 )   nstock = nlast + 1 
     377      ! 
     378   END SUBROUTINE sbc_lim_init 
     379 
     380 
     381   SUBROUTINE ice_run 
     382      !!------------------------------------------------------------------- 
     383      !!                  ***  ROUTINE ice_run *** 
     384      !!                  
     385      !! ** Purpose :   Definition some run parameter for ice model 
     386      !! 
     387      !! ** Method  :   Read the namicerun namelist and check the parameter  
     388      !!              values called at the first timestep (nit000) 
     389      !! 
     390      !! ** input   :   Namelist namicerun 
     391      !!------------------------------------------------------------------- 
     392      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, ln_nicep, ln_limdiahsb, ln_limdiaout 
     393      INTEGER  ::   ios                 ! Local integer output status for namelist read 
     394      !!------------------------------------------------------------------- 
     395      !                     
     396      REWIND( numnam_ice_ref )              ! Namelist namicerun in reference namelist : Parameters for ice 
     397      READ  ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901) 
     398901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp ) 
     399 
     400      REWIND( numnam_ice_cfg )              ! Namelist namicerun in configuration namelist : Parameters for ice 
     401      READ  ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 ) 
     402902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp ) 
     403      IF(lwm) WRITE ( numoni, namicerun ) 
     404      ! 
     405      ! 
     406      IF(lwp) THEN                        ! control print 
     407         WRITE(numout,*) 
     408         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 
     409         WRITE(numout,*) ' ~~~~~~' 
     410         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn 
     411         WRITE(numout,*) '   maximum ice concentration                               = ', amax  
     412         WRITE(numout,*) '   Several ice points in the ice or not in ocean.output    = ', ln_nicep 
     413         WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb 
     414         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout 
     415      ENDIF 
     416      ! 
     417      !IF( lk_mpp .AND. ln_nicep ) THEN 
     418      !   ln_nicep = .FALSE. 
     419      !   CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' ) 
     420      !ENDIF 
     421      IF( ln_nicep ) THEN      ! control print at a given point 
     422         jiindx = 15    ;   jjindx =  44 
     423         IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
     424      ENDIF 
     425      ! 
     426   END SUBROUTINE ice_run 
     427 
     428 
     429   SUBROUTINE lim_itd_init 
     430      !!------------------------------------------------------------------ 
     431      !!                ***  ROUTINE lim_itd_init *** 
     432      !! 
     433      !! ** Purpose :   Initializes the ice thickness distribution 
     434      !! ** Method  :   ... 
     435      !!    Note    : hi_max(jpl) is here set up to a value close to 7 m for 
     436      !!              limistate (only) and is changed to 99 m in sbc_lim_init 
     437      !!------------------------------------------------------------------ 
     438      INTEGER  ::   jl                   ! dummy loop index 
     439      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars 
     440      REAL(wp) ::   zhmax, znum, zden, zalpha ! 
     441      !!------------------------------------------------------------------ 
     442 
     443      IF(lwp) WRITE(numout,*) 
     444      IF(lwp) WRITE(numout,*) 'lim_itd_init : Initialization of ice thickness distribution ' 
     445      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     446 
     447      !------------------------------------------------------------------------------! 
     448      ! 1) Ice thickness distribution parameters initialization     
     449      !------------------------------------------------------------------------------! 
     450      IF(lwp) WRITE(numout,*) ' Number of ice categories jpl = ', jpl 
     451 
     452      !- Thickness categories boundaries  
     453      !---------------------------------- 
     454      !  Clem - je sais pas encore dans quelle namelist les mettre, ca depend des chgts liés à bdy 
     455      nn_hibnd  = 2                !  thickness category boundaries: tanh function (1) h^(-alpha) (2) 
     456      rn_hibnd  = 2.5              !  mean thickness of the domain (used to compute category boundaries, nn_hibnd = 2 only) 
     457 
     458      hi_max(:) = 0._wp 
     459 
     460      SELECT CASE ( nn_hibnd  )        
     461                                   !---------------------- 
     462         CASE (1)                  ! tanh function (CICE) 
     463                                   !---------------------- 
     464         zc1 =  3._wp / REAL( jpl, wp ) 
     465         zc2 = 10._wp * zc1 
     466         zc3 =  3._wp 
     467 
     468         DO jl = 1, jpl 
     469            zx1 = REAL( jl-1, wp ) / REAL( jpl, wp ) 
     470            hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) ) 
     471         END DO 
     472 
     473                                   !---------------------- 
     474         CASE (2)                  ! h^(-alpha) function 
     475                                   !---------------------- 
     476         zalpha = 0.05             ! exponent of the transform function 
     477 
     478         zhmax  = 3.*rn_hibnd 
     479 
     480         DO jl = 1, jpl  
     481            znum = jpl * ( zhmax+1 )**zalpha 
     482            zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl 
     483            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 
     484         END DO 
     485 
     486      END SELECT 
     487 
     488      DO jl = 1, jpl 
     489         hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 
     490      END DO 
     491      ! Set hi_max(jpl) to a big value to ensure that all ice is thinner than hi_max(jpl) 
     492      hi_max(jpl) = 99._wp 
     493 
     494      IF(lwp) WRITE(numout,*) ' Thickness category boundaries ' 
     495      IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl) 
     496      ! 
     497   END SUBROUTINE lim_itd_init 
     498 
    350499    
    351       SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   & 
     500   SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   & 
    352501         &                          pdqn_ice, pqla_ice, pdql_ice, k_limflx ) 
    353502      !!--------------------------------------------------------------------- 
    … …  
    427576      ! 
    428577   END SUBROUTINE ice_lim_flx 
    429     
    430     
    431    SUBROUTINE lim_ctl( kt ) 
    432       !!----------------------------------------------------------------------- 
    433       !!                   ***  ROUTINE lim_ctl ***  
    434       !!                  
    435       !! ** Purpose :   Alerts in case of model crash 
    436       !!------------------------------------------------------------------- 
    437       INTEGER, INTENT(in) ::   kt      ! ocean time step 
    438       INTEGER  ::   ji, jj, jk,  jl   ! dummy loop indices 
    439       INTEGER  ::   inb_altests       ! number of alert tests (max 20) 
    440       INTEGER  ::   ialert_id         ! number of the current alert 
    441       REAL(wp) ::   ztmelts           ! ice layer melting point 
    442       CHARACTER (len=30), DIMENSION(20)      ::   cl_alname   ! name of alert 
    443       INTEGER           , DIMENSION(20)      ::   inb_alp     ! number of alerts positive 
    444       !!------------------------------------------------------------------- 
    445  
    446       inb_altests = 10 
    447       inb_alp(:)  =  0 
    448  
    449       ! Alert if incompatible volume and concentration 
    450       ialert_id = 2 ! reference number of this alert 
    451       cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert 
    452  
    453       DO jl = 1, jpl 
    454          DO jj = 1, jpj 
    455             DO ji = 1, jpi 
    456                IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN 
    457                   !WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    458                   !WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
    459                   !WRITE(numout,*) ' Point - category', ji, jj, jl 
    460                   !WRITE(numout,*) ' a_i *** a_i_b   ', a_i      (ji,jj,jl), a_i_b  (ji,jj,jl) 
    461                   !WRITE(numout,*) ' v_i *** v_i_b   ', v_i      (ji,jj,jl), v_i_b  (ji,jj,jl) 
    462                   !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl) 
    463                   !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl) 
    464                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    465                ENDIF 
    466             END DO 
    467          END DO 
    468       END DO 
    469  
    470       ! Alerte if very thick ice 
    471       ialert_id = 3 ! reference number of this alert 
    472       cl_alname(ialert_id) = ' Very thick ice               ' ! name of the alert 
    473       jl = jpl  
    474       DO jj = 1, jpj 
    475          DO ji = 1, jpi 
    476             IF(   ht_i(ji,jj,jl)  >  50._wp   ) THEN 
    477                !CALL lim_prt_state( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
    478                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    479             ENDIF 
    480          END DO 
    481       END DO 
    482  
    483       ! Alert if very fast ice 
    484       ialert_id = 4 ! reference number of this alert 
    485       cl_alname(ialert_id) = ' Very fast ice               ' ! name of the alert 
    486       DO jj = 1, jpj 
    487          DO ji = 1, jpi 
    488             IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  & 
    489                &  at_i(ji,jj) > 0._wp   ) THEN 
    490                !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
    491                !WRITE(numout,*) ' ice strength             : ', strength(ji,jj) 
    492                !WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj)  
    493                !WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj) 
    494                !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)  
    495                !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj) 
    496                !WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj) 
    497                !WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj) 
    498                !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj) 
    499                !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj) 
    500                !WRITE(numout,*)  
    501                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    502             ENDIF 
    503          END DO 
    504       END DO 
    505  
    506       ! Alert if there is ice on continents 
    507       ialert_id = 6 ! reference number of this alert 
    508       cl_alname(ialert_id) = ' Ice on continents           ' ! name of the alert 
    509       DO jj = 1, jpj 
    510          DO ji = 1, jpi 
    511             IF(   tms(ji,jj) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN  
    512                !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
    513                !WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj)  
    514                !WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    515                !WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    516                !WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj) 
    517                !WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj) 
    518                !WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1) 
    519                !WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj) 
    520                !WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj) 
    521                ! 
    522                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    523             ENDIF 
    524          END DO 
    525       END DO 
    526  
    527 ! 
    528 !     ! Alert if very fresh ice 
    529       ialert_id = 7 ! reference number of this alert 
    530       cl_alname(ialert_id) = ' Very fresh ice               ' ! name of the alert 
    531       DO jl = 1, jpl 
    532          DO jj = 1, jpj 
    533             DO ji = 1, jpi 
    534                IF( sm_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    535 !                 CALL lim_prt_state(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
    536 !                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    537 !                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    538 !                 WRITE(numout,*)  
    539                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    540                ENDIF 
    541             END DO 
    542          END DO 
    543       END DO 
    544 ! 
    545  
    546 !     ! Alert if too old ice 
    547       ialert_id = 9 ! reference number of this alert 
    548       cl_alname(ialert_id) = ' Very old   ice               ' ! name of the alert 
    549       DO jl = 1, jpl 
    550          DO jj = 1, jpj 
    551             DO ji = 1, jpi 
    552                IF ( ( ( ABS( o_i(ji,jj,jl) ) > rdt_ice ) .OR. & 
    553                       ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 
    554                              ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
    555                   !CALL lim_prt_state( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
    556                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    557                ENDIF 
    558             END DO 
    559          END DO 
    560       END DO 
    561   
    562       ! Alert on salt flux 
    563       ialert_id = 5 ! reference number of this alert 
    564       cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert 
    565       DO jj = 1, jpj 
    566          DO ji = 1, jpi 
    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) , ' a_i_b      : ', a_i_b  (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) , ' v_i_b      : ', v_i_b  (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 
    577                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    578             ENDIF 
    579          END DO 
    580       END DO 
    581  
    582       ! Alert if qns very big 
    583       ialert_id = 8 ! reference number of this alert 
    584       cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert 
    585       DO jj = 1, jpj 
    586          DO ji = 1, jpi 
    587             IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
    588                ! 
    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                ! 
    595                !CALL lim_prt_state( kt, ji, jj, 2, '   ') 
    596                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    597                ! 
    598             ENDIF 
    599          END DO 
    600       END DO 
    601       !+++++ 
    602   
    603       ! Alert if very warm ice 
    604       ialert_id = 10 ! reference number of this alert 
    605       cl_alname(ialert_id) = ' Very warm ice                ' ! name of the alert 
    606       inb_alp(ialert_id) = 0 
    607       DO jl = 1, jpl 
    608          DO jk = 1, nlay_i 
    609             DO jj = 1, jpj 
    610                DO ji = 1, jpi 
    611                   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-10   & 
    613                      &                             .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 
    620                      inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    621                   ENDIF 
    622                END DO 
    623             END DO 
    624          END DO 
    625       END DO 
    626  
    627       ! sum of the alerts on all processors 
    628       IF( lk_mpp ) THEN 
    629          DO ialert_id = 1, inb_altests 
    630             CALL mpp_sum(inb_alp(ialert_id)) 
    631          END DO 
    632       ENDIF 
    633  
    634       ! print alerts 
    635       IF( lwp ) THEN 
    636          ialert_id = 1                                 ! reference number of this alert 
    637          cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert 
    638          WRITE(numout,*) ' time step ',kt 
    639          WRITE(numout,*) ' All alerts at the end of ice model ' 
    640          DO ialert_id = 1, inb_altests 
    641             WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! ' 
    642          END DO 
    643       ENDIF 
    644      ! 
    645    END SUBROUTINE lim_ctl 
    646   
    647     
    648    SUBROUTINE lim_prt_state( kt, ki, kj, kn, cd1 ) 
    649       !!----------------------------------------------------------------------- 
    650       !!                   ***  ROUTINE lim_prt_state ***  
    651       !!                  
    652       !! ** Purpose :   Writes global ice state on the (i,j) point  
    653       !!                in ocean.ouput  
    654       !!                3 possibilities exist  
    655       !!                n = 1/-1 -> simple ice state (plus Mechanical Check if -1) 
    656       !!                n = 2    -> exhaustive state 
    657       !!                n = 3    -> ice/ocean salt fluxes 
    658       !! 
    659       !! ** input   :   point coordinates (i,j)  
    660       !!                n : number of the option 
    661       !!------------------------------------------------------------------- 
    662       INTEGER         , INTENT(in) ::   kt            ! ocean time step 
    663       INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices 
    664       CHARACTER(len=*), INTENT(in) ::   cd1           ! 
    665       !! 
    666       INTEGER :: jl, ji, jj 
    667       !!------------------------------------------------------------------- 
    668  
    669       DO ji = mi0(ki), mi1(ki) 
    670          DO jj = mj0(kj), mj1(kj) 
    671  
    672             WRITE(numout,*) ' time step ',kt,' ',cd1             ! print title 
    673  
    674             !---------------- 
    675             !  Simple state 
    676             !---------------- 
    677              
    678             IF ( kn == 1 .OR. kn == -1 ) THEN 
    679                WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
    680                WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    681                WRITE(numout,*) ' Simple state ' 
    682                WRITE(numout,*) ' masks s,u,v   : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 
    683                WRITE(numout,*) ' lat - long    : ', gphit(ji,jj), glamt(ji,jj) 
    684                WRITE(numout,*) ' Time step     : ', numit 
    685                WRITE(numout,*) ' - Ice drift   ' 
    686                WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    687                WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj) 
    688                WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj) 
    689                WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1) 
    690                WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
    691                WRITE(numout,*) ' strength      : ', strength(ji,jj) 
    692                WRITE(numout,*) 
    693                WRITE(numout,*) ' - Cell values ' 
    694                WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    695                WRITE(numout,*) ' cell area     : ', area(ji,jj) 
    696                WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    697                WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
    698                WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
    699                DO jl = 1, jpl 
    700                   WRITE(numout,*) ' - Category (', jl,')' 
    701                   WRITE(numout,*) ' a_i           : ', a_i(ji,jj,jl) 
    702                   WRITE(numout,*) ' ht_i          : ', ht_i(ji,jj,jl) 
    703                   WRITE(numout,*) ' ht_s          : ', ht_s(ji,jj,jl) 
    704                   WRITE(numout,*) ' v_i           : ', v_i(ji,jj,jl) 
    705                   WRITE(numout,*) ' v_s           : ', v_s(ji,jj,jl) 
    706                   WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)/1.0e9 
    707                   WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)/1.0e9 
    708                   WRITE(numout,*) ' t_su          : ', t_su(ji,jj,jl) 
    709                   WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1,jl) 
    710                   WRITE(numout,*) ' t_i           : ', t_i(ji,jj,1:nlay_i,jl) 
    711                   WRITE(numout,*) ' sm_i          : ', sm_i(ji,jj,jl) 
    712                   WRITE(numout,*) ' smv_i         : ', smv_i(ji,jj,jl) 
    713                   WRITE(numout,*) 
    714                END DO 
    715             ENDIF 
    716             IF( kn == -1 ) THEN 
    717                WRITE(numout,*) ' Mechanical Check ************** ' 
    718                WRITE(numout,*) ' Check what means ice divergence ' 
    719                WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj) 
    720                WRITE(numout,*) ' Total lead fraction     ', ato_i(ji,jj) 
    721                WRITE(numout,*) ' Sum of both             ', ato_i(ji,jj) + at_i(ji,jj) 
    722                WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ji,jj) + at_i(ji,jj) - 1.00 
    723             ENDIF 
    724              
    725  
    726             !-------------------- 
    727             !  Exhaustive state 
    728             !-------------------- 
    729              
    730             IF ( kn .EQ. 2 ) THEN 
    731                WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
    732                WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    733                WRITE(numout,*) ' Exhaustive state ' 
    734                WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 
    735                WRITE(numout,*) ' Time step ', numit 
    736                WRITE(numout,*)  
    737                WRITE(numout,*) ' - Cell values ' 
    738                WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    739                WRITE(numout,*) ' cell area     : ', area(ji,jj) 
    740                WRITE(numout,*) ' at_i          : ', at_i(ji,jj)        
    741                WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)        
    742                WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)        
    743                WRITE(numout,*) ' u_ice(i-1,j)  : ', u_ice(ji-1,jj) 
    744                WRITE(numout,*) ' u_ice(i  ,j)  : ', u_ice(ji,jj) 
    745                WRITE(numout,*) ' v_ice(i  ,j-1): ', v_ice(ji,jj-1) 
    746                WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
    747                WRITE(numout,*) ' strength      : ', strength(ji,jj) 
    748                WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj) 
    749                WRITE(numout,*) ' u_ice_b       : ', u_ice_b(ji,jj)    , ' v_ice_b       : ', v_ice_b(ji,jj)   
    750                WRITE(numout,*) 
    751                 
    752                DO jl = 1, jpl 
    753                   WRITE(numout,*) ' - Category (',jl,')' 
    754                   WRITE(numout,*) '   ~~~~~~~~         '  
    755                   WRITE(numout,*) ' ht_i       : ', ht_i(ji,jj,jl)             , ' ht_s       : ', ht_s(ji,jj,jl) 
    756                   WRITE(numout,*) ' t_i        : ', t_i(ji,jj,1:nlay_i,jl) 
    757                   WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl) 
    758                   WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl) 
    759                   WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' a_i_b      : ', a_i_b(ji,jj,jl)    
    760                   WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    761                   WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' v_i_b      : ', v_i_b(ji,jj,jl)    
    762                   WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    763                   WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' v_s_b      : ', v_s_b(ji,jj,jl)   
    764                   WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl) 
    765                   WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' ei1        : ', e_i_b(ji,jj,1,jl)/1.0e9  
    766                   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 
    767                   WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' ei2_b      : ', e_i_b(ji,jj,2,jl)/1.0e9   
    768                   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 
    769                   WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl)  
    770                   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) 
    771                   WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' smv_i_b    : ', smv_i_b(ji,jj,jl)    
    772                   WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl)  
    773                   WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' oa_i_b     : ', oa_i_b(ji,jj,jl) 
    774                   WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl) 
    775                END DO !jl 
    776                 
    777                WRITE(numout,*) 
    778                WRITE(numout,*) ' - Heat / FW fluxes ' 
    779                WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    780                WRITE(numout,*) ' - Heat fluxes in and out the ice ***' 
    781                WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) ) 
    782                WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) ) 
    783                WRITE(numout,*) 
    784                WRITE(numout,*)  
    785                WRITE(numout,*) ' sst        : ', sst_m(ji,jj)   
    786                WRITE(numout,*) ' sss        : ', sss_m(ji,jj)   
    787                WRITE(numout,*)  
    788                WRITE(numout,*) ' - Stresses ' 
    789                WRITE(numout,*) '   ~~~~~~~~ ' 
    790                WRITE(numout,*) ' utau_ice   : ', utau_ice(ji,jj)  
    791                WRITE(numout,*) ' vtau_ice   : ', vtau_ice(ji,jj) 
    792                WRITE(numout,*) ' utau       : ', utau    (ji,jj)  
    793                WRITE(numout,*) ' vtau       : ', vtau    (ji,jj) 
    794                WRITE(numout,*) ' oc. vel. u : ', u_oce   (ji,jj) 
    795                WRITE(numout,*) ' oc. vel. v : ', v_oce   (ji,jj) 
    796             ENDIF 
    797              
    798             !--------------------- 
    799             ! Salt / heat fluxes 
    800             !--------------------- 
    801              
    802             IF ( kn .EQ. 3 ) THEN 
    803                WRITE(numout,*) ' lim_prt_state - Point : ',ji,jj 
    804                WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    805                WRITE(numout,*) ' - Salt / Heat Fluxes ' 
    806                WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    807                WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj) 
    808                WRITE(numout,*) ' Time step ', numit 
    809                WRITE(numout,*) 
    810                WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
    811                WRITE(numout,*) ' qsr       : ', qsr(ji,jj) 
    812                WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    813                WRITE(numout,*) 
    814                WRITE(numout,*) ' hfx_mass     : ', hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_snw(ji,jj) + hfx_res(ji,jj) 
    815                WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj) 
    816                WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj) 
    817                WRITE(numout,*) ' dhc          : ', diag_heat_dhc(ji,jj)               
    818                WRITE(numout,*) 
    819                WRITE(numout,*) ' hfx_dyn      : ', hfx_dyn(ji,jj) 
    820                WRITE(numout,*) ' hfx_thd      : ', hfx_thd(ji,jj) 
    821                WRITE(numout,*) ' hfx_res      : ', hfx_res(ji,jj) 
    822                WRITE(numout,*) ' fhtur        : ', fhtur(ji,jj)  
    823                WRITE(numout,*) ' qlead        : ', qlead(ji,jj) * r1_rdtice 
    824                WRITE(numout,*) 
    825                WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
    826                WRITE(numout,*) ' emp       : ', emp    (ji,jj) 
    827                WRITE(numout,*) ' sfx       : ', sfx    (ji,jj) 
    828                WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj) 
    829                WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
    830                WRITE(numout,*) ' sfx_dyn   : ', sfx_dyn(ji,jj) 
    831                WRITE(numout,*) 
    832                WRITE(numout,*) ' - Momentum fluxes ' 
    833                WRITE(numout,*) ' utau      : ', utau(ji,jj)  
    834                WRITE(numout,*) ' vtau      : ', vtau(ji,jj) 
    835             ENDIF  
    836             WRITE(numout,*) ' ' 
    837             ! 
    838          END DO 
    839       END DO 
    840       ! 
    841    END SUBROUTINE lim_prt_state 
    842     
    843       
     578 
     579   SUBROUTINE sbc_lim_flx0 
     580      !!---------------------------------------------------------------------- 
     581      !!                  ***  ROUTINE sbc_lim_flx0  *** 
     582      !! 
     583      !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining 
     584      !!               of the time step 
     585      !!---------------------------------------------------------------------- 
     586      sfx    (:,:) = 0._wp   ; 
     587      sfx_bri(:,:) = 0._wp   ;  
     588      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
     589      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
     590      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
     591      sfx_res(:,:) = 0._wp 
     592       
     593      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
     594      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
     595      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
     596      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
     597      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
     598      wfx_spr(:,:) = 0._wp   ;    
     599       
     600      hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
     601      hfx_thd(:,:) = 0._wp   ;    
     602      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
     603      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
     604      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
     605      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
     606      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
     607      hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
     608 
     609      afx_tot(:,:) = 0._wp   ; 
     610      afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp 
     611       
     612   END SUBROUTINE sbc_lim_flx0 
     613       
    844614   FUNCTION fice_cell_ave ( ptab ) 
    845615      !!-------------------------------------------------------------------------- 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4990 r5051  
    271271      IF( ln_ssr           )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation 
    272272      ! 
     273      IF( nn_ice == 3      )   CALL sbc_lim_init               ! LIM3 initialisation 
     274 
    273275      IF( nn_ice == 4      )   CALL cice_sbc_init( nsbc )      ! CICE initialisation 
    274276      ! 
Note: See TracChangeset for help on using the changeset viewer.