Changeset 8522


Ignore:
Timestamp:
2017-09-14T17:52:02+02:00 (3 years ago)
Author:
clem
Message:

changes in style - part6 - ice diffusion (hope the scheme still converges)

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90

    r8518 r8522  
    1919   USE phycst         ! physical constant 
    2020   USE daymod         ! model calendar 
    21    USE sbc_oce , ONLY : sfx   ! surface boundary condition: ocean fields 
     21   USE sbc_oce , ONLY : sfx, nn_fsbc   ! surface boundary condition: ocean fields 
    2222   USE icerst         ! ice restart 
    2323   ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

    r8518 r8522  
    280280      !!------------------------------------------------------------------- 
    281281      INTEGER  ::   ji, jk   ! dummy loop indices 
    282       REAL(wp) ::   ztmelts, z1_2cp, zbbb, zccc  ! local scalar  
     282      REAL(wp) ::   ztmelts, zbbb, zccc  ! local scalar  
    283283      !!------------------------------------------------------------------- 
    284284      ! Recover ice temperature 
    285       z1_2cp = 1._wp / ( 2._wp * cpic ) 
    286285      DO jk = 1, nlay_i 
    287286         DO ji = 1, nidx 
     
    290289            zbbb          = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus 
    291290            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) ) 
    292             t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * z1_2cp 
     291            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpic 
    293292             
    294293            ! mask temperature 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

    r8514 r8522  
    7878      REAL(wp) ::   zgrr         ! bottom growth rate 
    7979      REAL(wp) ::   zt_i_new     ! bottom formation temperature 
     80      REAL(wp) ::   z1_rho       ! 1/(rhosn+rau0-rhoic) 
    8081 
    8182      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean 
     
    596597      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,  
    597598      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 
     599      z1_rho = 1._wp / ( rhosn+rau0-rhoic ) 
    598600      DO ji = 1, nidx 
    599601         ! 
    600          dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) 
     602         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) * z1_rho ) 
    601603 
    602604         ht_i_1d(ji)    = ht_i_1d(ji) + dh_snowice(ji) 
     
    646648            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 
    647649            ! recalculate t_s_1d from e_s_1d 
    648             t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 
     650            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhosn * r1_cpic + lfus * r1_cpic ) 
    649651         END DO 
    650652      END DO 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90

    r8518 r8522  
    7272      !!           ice salinities          : s_i_1d 
    7373      !!           number of layers in the ice/snow: nlay_i, nlay_s 
    74       !!           profile of the ice/snow layers : z_i, z_s 
    7574      !!           total ice/snow thickness : ht_i_1d, ht_s_1d 
    7675      !!------------------------------------------------------------------ 
    77       INTEGER ::   ji             ! spatial loop index 
    78       INTEGER ::   ii, ij         ! temporary dummy loop index 
     76      INTEGER ::   ji, jk         ! spatial loop index 
    7977      INTEGER ::   numeq          ! current reference number of equation 
    80       INTEGER ::   jk             ! vertical dummy loop index  
    8178      INTEGER ::   minnumeqmin, maxnumeqmax 
    8279      INTEGER ::   iconv          ! number of iterations in iterative procedure 
     
    9390      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity 
    9491      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C  
     92      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature  
    9593      REAL(wp) ::   ztmelt_i                  ! ice melting temperature 
    9694      REAL(wp) ::   z1_hsu 
    9795      REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
    98       REAL(wp) ::   zdti_bnd = 1.e-4_wp       ! maximal authorized error on temperature  
    9996      REAL(wp) ::   zcpi                      ! Ice specific heat 
    100       REAL(wp) ::   zhi                       !  
    10197      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat 
     98      REAL(wp) ::   zfac                      ! dummy factor 
    10299       
    103100      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow 
    104101      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration 
    105       REAL(wp), DIMENSION(jpij)     ::   zh_i        ! ice layer thickness 
    106       REAL(wp), DIMENSION(jpij)     ::   zh_s        ! snow layer thickness 
     102      REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness 
     103      REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness 
    107104      REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface 
    108105      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface 
    109106      REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function 
    110107      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function 
    111       REAL(wp), DIMENSION(jpij)     ::   zdti        ! current error on temperature 
    112108      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice 
    113109       
    114       REAL(wp), DIMENSION(jpij,nlay_i)   ::   z_i         ! Vertical cotes of the layers in the ice 
    115       REAL(wp), DIMENSION(jpij,nlay_s)   ::   z_s         ! Vertical cotes of the layers in the snow 
    116       REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztib        ! Old temperature in the ice 
    117       REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsb        ! Temporary temperature in the snow 
    118       REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
    119       REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
     110      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice 
     111      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow 
     112      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence 
     113      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence 
    120114      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity 
    121115      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
     
    136130      ! Mono-category 
    137131      REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done 
    138       REAL(wp) :: zratio_s      ! dummy factor 
    139       REAL(wp) :: zratio_i      ! dummy factor 
    140       REAL(wp) :: zh_thres      ! thickness thres. for G(h) computation 
    141132      REAL(wp) :: zhe           ! dummy factor 
    142       REAL(wp) :: zkimean       ! mean sea ice thermal conductivity 
    143       REAL(wp) :: zfac          ! dummy factor 
    144       REAL(wp) :: zihe          ! dummy factor 
    145       REAL(wp) :: zheshth       ! dummy factor 
     133      REAL(wp) :: zcnd_i        ! mean sea ice thermal conductivity 
    146134      !!------------------------------------------------------------------      
    147135 
     
    155143      ! 1) Initialization                                                            ! 
    156144      !------------------------------------------------------------------------------! 
    157       DO ji = 1 , nidx 
     145      DO ji = 1, nidx 
    158146         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not 
    159147         ! layer thickness 
     
    161149         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 
    162150      END DO 
    163  
    164       !-------------------- 
    165       ! Ice / snow layers 
    166       !-------------------- 
    167       z_s(1:nidx,1) = zh_s(1:nidx) 
    168       z_i(1:nidx,1) = zh_i(1:nidx) 
    169       DO jk = 2, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer 
    170          DO ji = 1 , nidx 
    171             z_s(ji,jk) = z_s(ji,jk-1) + zh_s(ji) 
    172          END DO 
    173       END DO 
    174  
    175       DO jk = 2, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer 
    176          DO ji = 1 , nidx 
    177             z_i(ji,jk) = z_i(ji,jk-1) + zh_i(ji) 
    178          END DO 
    179       END DO 
     151      ! 
     152      WHERE( zh_i(1:nidx) >= epsi10 )   ;   z1_h_i(1:nidx) = 1._wp / zh_i(1:nidx) 
     153      ELSEWHERE                         ;   z1_h_i(1:nidx) = 0._wp 
     154      END WHERE 
     155 
     156      WHERE( zh_s(1:nidx) >= epsi10 )   ;   z1_h_s(1:nidx) = 1._wp / zh_s(1:nidx) 
     157      ELSEWHERE                         ;   z1_h_s(1:nidx) = 0._wp 
     158      END WHERE 
     159      ! 
     160      ! temperatures 
     161      ztsub  (1:nidx)   = t_su_1d(1:nidx)  ! temperature at the previous iteration 
     162      ztsold (1:nidx,:) = t_s_1d(1:nidx,:)   ! Old snow temperature 
     163      ztiold (1:nidx,:) = t_i_1d(1:nidx,:)   ! Old ice temperature 
     164      t_su_1d(1:nidx)   = MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary 
    180165      ! 
    181166      !------------------------------------------------------------------------------| 
    182       ! 2) Radiation                                                       | 
     167      ! 2) Radiation                                                              | 
    183168      !------------------------------------------------------------------------------| 
    184169      ! 
    185170      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 
    186       DO ji = 1 , nidx 
     171      DO ji = 1, nidx 
    187172         !------------------- 
    188173         ! Computation of i0 
     
    196181         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
    197182         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 
    198          zhi = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )      
    199          i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zhi * fr2_i0_1d(ji) ) 
     183         zfac = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )      
     184         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) ) 
    200185 
    201186         !------------------------------------------------------- 
     
    216201         DO ji = 1, nidx 
    217202            !                             ! radiation transmitted below the layer-th snow layer 
    218             zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * z_s(ji,jk) ) 
     203            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) ) 
    219204            !                             ! radiation absorbed by the layer-th snow layer 
    220205            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 
     
    226211         DO ji = 1, nidx 
    227212            !                             ! radiation transmitted below the layer-th ice layer 
    228             zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * z_i(ji,jk) ) 
     213            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) 
    229214            !                             ! radiation absorbed by the layer-th ice layer 
    230215            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 
     
    238223      !------------------------------------------------------------------------------| 
    239224      ! 
    240       ztsub  (1:nidx)   =  t_su_1d(1:nidx)                          ! temperature at the previous iter 
    241       t_su_1d(1:nidx)   =  MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary 
    242       zdti   (1:nidx)   =  1000._wp                                 ! initial value of error 
    243       ztsb   (1:nidx,:) =  t_s_1d(1:nidx,:)                         ! Old snow temperature 
    244       ztib   (1:nidx,:) =  t_i_1d(1:nidx,:)                         ! Old ice temperature 
    245  
    246       iconv    =  0           ! number of iterations 
    247       zdti_max =  1000._wp    ! maximal value of error on all points 
    248  
     225      iconv    =  0          ! number of iterations 
     226      zdti_max =  1000._wp   ! maximal value of error on all points 
    249227      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) 
    250228         ! 
    251229         iconv = iconv + 1 
    252230         ! 
     231         ztib(1:nidx,:) = t_i_1d(1:nidx,:) 
     232         ztsb(1:nidx,:) = t_s_1d(1:nidx,:) 
     233         ! 
    253234         !------------------------------------------------------------------------------| 
    254235         ! 4) Sea ice thermal conductivity                                              | 
     
    256237         ! 
    257238         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T 
    258             DO ji = 1 , nidx 
    259                ztcond_i(ji,0) = MAX( zkimin, rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) ) 
     239            ! 
     240            DO ji = 1, nidx 
     241               ztcond_i(ji,0)      = rcdic + zbeta * s_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
     242               ztcond_i(ji,nlay_i) = rcdic + zbeta * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
    260243            END DO 
    261244            DO jk = 1, nlay_i-1 
    262                DO ji = 1 , nidx 
    263                   ztcond_i(ji,jk) = MAX( zkimin, rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  & 
    264                      &                           MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) ) 
     245               DO ji = 1, nidx 
     246                  ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  & 
     247                     &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
    265248               END DO 
    266249            END DO 
    267  
     250            ! 
    268251         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 
    269             DO ji = 1 , nidx 
    270                ztcond_i(ji,0) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   & 
    271                   &                               - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) ) 
     252            ! 
     253            DO ji = 1, nidx 
     254               ztcond_i(ji,0)      = rcdic + 0.09_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
     255                  &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
     256               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )  & 
     257                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
    272258            END DO 
    273259            DO jk = 1, nlay_i-1 
    274                DO ji = 1 , nidx 
    275                   ztcond_i(ji,jk) = MAX( zkimin, rcdic +                                                           &  
    276                      &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              & 
    277                      &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   & 
    278                      &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )  ) 
     260               DO ji = 1, nidx 
     261                  ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /          & 
     262                     &                     MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )   & 
     263                     &                    - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
    279264               END DO 
    280265            END DO 
    281             DO ji = 1 , nidx 
    282                ztcond_i(ji,nlay_i) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   & 
    283                   &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) )   
    284             END DO 
     266            ! 
    285267         ENDIF 
    286           
     268         ztcond_i(1:nidx,:) = MAX( zkimin, ztcond_i(1:nidx,:) )         
    287269         ! 
    288270         !------------------------------------------------------------------------------| 
     
    294276         ! Fichefet and Morales Maqueda, JGR 1997 
    295277 
    296          zghe(:) = 1._wp 
    297  
     278         zghe(1:nidx) = 1._wp 
     279          
    298280         SELECT CASE ( nn_monocat ) 
    299281 
    300          CASE (1,3) 
     282         CASE ( 1 , 3 ) 
    301283 
    302284            zepsilon = 0.1_wp 
    303             zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp 
    304  
    305285            DO ji = 1, nidx 
    306286    
    307287               ! Mean sea ice thermal conductivity 
    308                zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp ) 
     288               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) 
    309289 
    310290               ! Effective thickness he (zhe) 
    311                zfac     = 1._wp / ( rn_cnd_s + zkimean ) 
    312                zratio_s = rn_cnd_s   * zfac 
    313                zratio_i = zkimean * zfac 
    314                zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 
     291               zhe  = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) 
    315292 
    316293               ! G(he) 
    317                rswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if > 
    318                zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) 
     294               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN 
     295                  zghe(ji) =  MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) ) 
     296               ENDIF 
    319297    
    320                ! Impose G(he) < 2. 
    321                zghe(ji) = MIN( zghe(ji), 2._wp ) 
    322  
    323298            END DO 
    324299 
    325300         END SELECT 
    326  
    327301         ! 
    328302         !------------------------------------------------------------------------------| 
     
    331305         ! 
    332306         !--- Snow 
    333          DO ji = 1, nidx 
    334             zfac                  =  1. / MAX( epsi10 , zh_s(ji) ) 
    335             zkappa_s(ji,0)        = zghe(ji) * rn_cnd_s * zfac 
    336             zkappa_s(ji,nlay_s)   = zghe(ji) * rn_cnd_s * zfac 
    337          END DO 
    338  
    339          DO jk = 1, nlay_s-1 
    340             DO ji = 1 , nidx 
    341                zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rn_cnd_s / MAX( epsi10, 2.0 * zh_s(ji) ) 
    342             END DO 
     307         DO jk = 0, nlay_s-1 
     308            DO ji = 1, nidx 
     309               zkappa_s(ji,jk)  = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
     310            END DO 
     311         END DO 
     312         DO ji = 1, nidx   ! Snow-ice interface 
     313            zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 
     314            IF( zfac > epsi10 ) THEN 
     315               zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 
     316            ELSE 
     317               zkappa_s(ji,nlay_s) = 0._wp 
     318            ENDIF 
    343319         END DO 
    344320 
    345321         !--- Ice 
    346          DO jk = 1, nlay_i-1 
    347             DO ji = 1 , nidx 
    348                zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) ) 
    349             END DO 
    350          END DO 
    351  
    352          !--- Snow-ice interface 
    353          DO ji = 1 , nidx 
    354             zfac                  = 1./ MAX( epsi10 , zh_i(ji) ) 
    355             zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac 
    356             zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 
    357             zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rn_cnd_s * ztcond_i(ji,0) / &  
    358            &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cnd_s * zh_i(ji) ) ) 
    359             zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
    360          END DO 
    361  
     322         DO jk = 0, nlay_i 
     323            DO ji = 1, nidx 
     324               zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 
     325            END DO 
     326         END DO 
     327         DO ji = 1, nidx   ! Snow-ice interface 
     328            zkappa_i(ji,0)     = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
     329         END DO 
    362330         ! 
    363331         !------------------------------------------------------------------------------| 
     
    366334         ! 
    367335         DO jk = 1, nlay_i 
    368             DO ji = 1 , nidx 
    369                ztitemp(ji,jk) = t_i_1d(ji,jk) 
    370                zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 
    371                zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zcpi * zh_i(ji), epsi10 ) 
     336            DO ji = 1, nidx 
     337               zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 
     338               zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi )  
    372339            END DO 
    373340         END DO 
    374341 
    375342         DO jk = 1, nlay_s 
    376             DO ji = 1 , nidx 
    377                ztstemp(ji,jk) = t_s_1d(ji,jk) 
    378                zeta_s(ji,jk)  = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 ) 
    379             END DO 
    380          END DO 
    381  
     343            DO ji = 1, nidx 
     344               zeta_s(ji,jk)  = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 
     345            END DO 
     346         END DO 
    382347         ! 
    383348         !------------------------------------------------------------------------------| 
     
    386351         ! 
    387352         IF ( ln_dqns_i ) THEN  
    388             DO ji = 1 , nidx 
     353            DO ji = 1, nidx 
    389354               ! update of the non solar flux according to the update in T_su 
    390355               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
     
    392357         ENDIF 
    393358 
    394          ! Update incoming flux 
    395          DO ji = 1 , nidx 
    396             ! update incoming flux 
    397             zf(ji)    =          zfsw(ji)              & ! net absorbed solar radiation 
    398                &         + qns_ice_1d(ji)                ! non solar total flux (LWup, LWdw, SH, LH) 
    399          END DO 
    400  
     359         DO ji = 1, nidx 
     360            zf(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 
     361         END DO 
    401362         ! 
    402363         !------------------------------------------------------------------------------| 
     
    412373 
    413374         DO numeq=1,nlay_i+3 
    414             DO ji = 1 , nidx 
     375            DO ji = 1, nidx 
    415376               ztrid(ji,numeq,1) = 0. 
    416377               ztrid(ji,numeq,2) = 0. 
     
    423384 
    424385         DO numeq = nlay_s + 2, nlay_s + nlay_i  
    425             DO ji = 1 , nidx 
     386            DO ji = 1, nidx 
    426387               jk                 = numeq - nlay_s - 1 
    427388               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 
    428389               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
    429390               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk) 
    430                zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
     391               zindterm(ji,numeq) =  ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
    431392            END DO 
    432393         ENDDO 
    433394 
    434395         numeq =  nlay_s + nlay_i + 1 
    435          DO ji = 1 , nidx 
     396         DO ji = 1, nidx 
    436397            !!ice bottom term 
    437398            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
    438399            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 
    439400            ztrid(ji,numeq,3)  =  0.0 
    440             zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
     401            zindterm(ji,numeq) =  ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
    441402               &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
    442403         ENDDO 
    443404 
    444405 
    445          DO ji = 1 , nidx 
     406         DO ji = 1, nidx 
    446407            IF ( ht_s_1d(ji) > 0.0 ) THEN 
    447408               ! 
     
    456417                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
    457418                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk) 
    458                   zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
     419                  zindterm(ji,numeq)  =  ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
    459420               END DO 
    460421 
    461422               !!case of only one layer in the ice (ice equation is altered) 
    462                IF ( nlay_i.eq.1 ) THEN 
     423               IF ( nlay_i == 1 ) THEN 
    463424                  ztrid(ji,nlay_s+2,3)    =  0.0 
    464425                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
     
    483444                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    484445                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1) 
    485                   zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
     446                  zindterm(ji,2) =  ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
    486447 
    487448               ELSE  
     
    498459                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    499460                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1)  
    500                   zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *  & 
     461                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  & 
    501462                     &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
    502463               ENDIF 
     
    526487                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
    527488                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1)   
    528                   zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)   
     489                  zindterm(ji,numeqmin(ji)+1)=  ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)   
    529490 
    530491                  !!case of only one layer in the ice (surface & ice equations are altered) 
     
    538499                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0 
    539500 
    540                      zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1) *  & 
     501                     zindterm(ji,numeqmin(ji)+1) =  ztiold(ji,1) + zeta_i(ji,1) *  & 
    541502                        &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 
    542503                  ENDIF 
     
    556517                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
    557518                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1) 
    558                   zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1) *  & 
     519                  zindterm(ji,numeqmin(ji))     =  ztiold(ji,1) + zeta_i(ji,1) *  & 
    559520                     &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
    560521 
     
    564525                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    565526                     ztrid(ji,numeqmin(ji),3)  =  0.0 
    566                      zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  & 
     527                     zindterm(ji,numeqmin(ji)) =  ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  & 
    567528                        &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 
    568529                  ENDIF 
     
    572533 
    573534         END DO 
    574  
    575535         ! 
    576536         !------------------------------------------------------------------------------| 
     
    578538         !------------------------------------------------------------------------------| 
    579539         ! 
    580  
    581540         ! Solve the tridiagonal system with Gauss elimination method. 
    582541         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,  
     
    586545         minnumeqmin = nlay_i+5 
    587546 
    588          DO ji = 1 , nidx 
     547         DO ji = 1, nidx 
    589548            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji)) 
    590549            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2) 
     
    594553 
    595554         DO jk = minnumeqmin+1, maxnumeqmax 
    596             DO ji = 1 , nidx 
     555            DO ji = 1, nidx 
    597556               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 
    598557               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1) 
     
    601560         END DO 
    602561 
    603          DO ji = 1 , nidx 
     562         DO ji = 1, nidx 
    604563            ! ice temperatures 
    605             t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji)) 
     564            t_i_1d(ji,nlay_i) =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji)) 
    606565         END DO 
    607566 
    608567         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1 
    609             DO ji = 1 , nidx 
     568            DO ji = 1, nidx 
    610569               jk    =  numeq - nlay_s - 1 
    611570               t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq) 
     
    613572         END DO 
    614573 
    615          DO ji = 1 , nidx 
     574         DO ji = 1, nidx 
    616575            ! snow temperatures       
    617             IF (ht_s_1d(ji) > 0._wp) & 
    618                t_s_1d(ji,nlay_s)     =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  & 
    619                &                        / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) )  
    620  
     576            IF( ht_s_1d(ji) > 0._wp ) THEN 
     577               t_s_1d(ji,nlay_s) =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  & 
     578                  &                 / zdiagbis(ji,nlay_s+1) 
     579            ENDIF 
    621580            ! surface temperature 
    622             isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) ) 
    623581            ztsub(ji) = t_su_1d(ji) 
    624             IF( t_su_1d(ji) < rt0 ) & 
     582            IF( t_su_1d(ji) < rt0 ) THEN 
    625583               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  & 
    626                &             ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
     584                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 
     585            ENDIF 
    627586         END DO 
    628587         ! 
     
    632591         ! 
    633592         ! check that nowhere it has started to melt 
    634          ! zdti(ji) is a measure of error, it has to be under zdti_bnd 
    635          DO ji = 1 , nidx 
    636             t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  ) 
    637             zdti   (ji) =  ABS( t_su_1d(ji) - ztsub(ji) )      
     593         ! zdti_max is a measure of error, it has to be under zdti_bnd 
     594         zdti_max = 0._wp 
     595         DO ji = 1, nidx 
     596            t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
     597            zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )      
    638598         END DO 
    639599 
    640600         DO jk  =  1, nlay_s 
    641             DO ji = 1 , nidx 
    642                t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rt0 ), 190._wp ) 
    643                zdti  (ji)    = MAX( zdti(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) ) 
     601            DO ji = 1, nidx 
     602               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
     603               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
    644604            END DO 
    645605         END DO 
    646606 
    647607         DO jk  =  1, nlay_i 
    648             DO ji = 1 , nidx 
     608            DO ji = 1, nidx 
    649609               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0  
    650                t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp ) 
    651                zdti  (ji)    =  MAX( zdti(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) ) 
     610               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 
     611               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    652612            END DO 
    653613         END DO 
     
    655615         ! Compute spatial maximum over all errors 
    656616         ! note that this could be optimized substantially by iterating only the non-converging points 
    657          zdti_max = 0._wp 
    658          DO ji = 1, nidx 
    659             zdti_max = MAX( zdti_max, zdti(ji) )    
    660          END DO 
    661617         IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
    662618 
    663619      END DO  ! End of the do while iterative procedure 
    664  
    665       ! MV SIMIP 2016 
    666       !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    667       DO ji = 1, nidx 
    668          zfac        = 1. / MAX( epsi10 , rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) ) 
    669          IF( zh_s(ji) >= 1.e-3 ) THEN 
    670             t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) + & 
    671                &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac 
    672          ELSE 
    673             t_si_1d(ji) = t_su_1d(ji) 
    674          ENDIF 
    675       END DO 
    676       ! END MV SIMIP 2016 
    677620 
    678621      IF( ln_icectl .AND. lwp ) THEN 
     
    687630      DO ji = 1, nidx 
    688631         !                                ! surface ice conduction flux 
    689          isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) 
    690632         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
    691633            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
     
    693635         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 
    694636      END DO 
    695  
    696       ! MV SIMIP 2016 
    697       !--- Conduction fluxes (positive downwards) 
    698       DO ji = 1, nidx 
    699          diag_fc_bo_1d(ji)   = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 
    700          diag_fc_su_1d(ji)   = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji) 
    701       END DO 
    702       ! END MV SIMIP 2016 
    703637 
    704638      ! --- computes sea ice energy of melting compulsory for icethd_dh --- ! 
     
    730664 
    731665      END DO    
     666 
     667      ! --- Diagnostics SIMIP --- ! 
     668      DO ji = 1, nidx 
     669         !--- Conduction fluxes (positive downwards) 
     670         diag_fc_bo_1d(ji)   = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 
     671         diag_fc_su_1d(ji)   = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji) 
     672 
     673         !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
     674         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 
     675         IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 
     676            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) + & 
     677               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 
     678         ELSE 
     679            t_si_1d(ji) = t_su_1d(ji) 
     680         ENDIF 
     681      END DO 
    732682      ! 
    733683   END SUBROUTINE ice_thd_dif 
     
    748698      DO jk = 1, nlay_i             ! Sea ice energy of melting 
    749699         DO ji = 1, nidx 
    750             ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0 
    751             t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point 
    752                                                           !   (sometimes dif scheme produces abnormally high temperatures)    
    753             e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                           & 
    754                &                    + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) )   & 
    755                &                    - rcp  *         ( ztmelts-rt0 )  )  
     700            ztmelts      = - tmut  * s_i_1d(ji,jk) 
     701            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point 
     702                                                                !   (sometimes dif scheme produces abnormally high temperatures)    
     703            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           & 
     704               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   & 
     705               &                    - rcp  *   ztmelts ) 
    756706         END DO 
    757707      END DO 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

    r8517 r8522  
    150150      !!------------------------------------------------------------------ 
    151151      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    152       REAL(wp) ::   ze_i, z1_cp, z1_2cp             ! local scalars 
     152      REAL(wp) ::   ze_i             ! local scalars 
    153153      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      - 
    154154      REAL(wp) ::   zhmax, z1_zhmax, zsm_i          !   -      - 
     
    193193 
    194194      !------------------- 
    195       ! Ice temperature   [K]   (with a minimum value (rt0 - 100.) imposed everywhere) 
     195      ! Ice temperature   [K]   (with a minimum value (rt0 - 100.)) 
    196196      !------------------- 
    197197      zlay_i   = REAL( nlay_i , wp )    ! number of layers 
    198       z1_2cp  = 1._wp / ( 2._wp * cpic ) 
    199198      DO jl = 1, jpl 
    200199         DO jk = 1, nlay_i 
     
    208207                     zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus 
    209208                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) ) 
    210                      t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * z1_2cp , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts 
     209                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpic , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts 
    211210                     ! 
    212211                  ELSE                                !--- no ice 
    213                      t_i(ji,jj,jk,jl) = rt0 - 100._wp                                   ! impose 173.15 K (i.e. -100 C) 
    214 !!clem: I think we should impose rt0 instead 
     212                     t_i(ji,jj,jk,jl) = rt0 
    215213                  ENDIF 
    216214               END DO 
     
    220218 
    221219      !-------------------- 
    222       ! Snow temperature   [K]   (with a minimum value (rt0 - 100.) imposed everywhere) 
     220      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.)) 
    223221      !-------------------- 
    224222      zlay_s = REAL( nlay_s , wp ) 
    225       z1_cp  = 1._wp / cpic 
    226223      DO jk = 1, nlay_s 
    227224         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area 
    228             t_s(:,:,jk,:) = MAX( -100._wp , MIN( z1_cp * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0 
     225            t_s(:,:,jk,:) = MAX( -100._wp , MIN( r1_cpic * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0 
    229226         ELSEWHERE                           !--- no ice 
    230             t_s(:,:,jk,:) =  rt0 - 100._wp         ! impose 173.15 K (i.e. -100 C) 
     227            t_s(:,:,jk,:) = rt0 
    231228         END WHERE 
    232229      END DO 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90

    r8518 r8522  
    179179      ENDIF 
    180180 
    181       ztmp = rday / rhosn 
     181      ztmp = rday * r1_rhosn 
    182182      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow) 
    183183      CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt  
     
    237237      ! Add-ons for SIMIP 
    238238      !-------------------------------- 
    239       zrho1 = ( rau0 - rhoic ) / rau0; zrho2 = rhosn / rau0 
     239      zrho1 = ( rau0 - rhoic ) * r1_rau0; zrho2 = rhosn * r1_rau0 
    240240 
    241241      IF ( iom_use( "icepres"  ) ) CALL iom_put( "icepres"     , zswi(:,:)                     )                                ! Ice presence (1 or 0)  
     
    278278 
    279279      IF ( iom_use( "dmsspr"   ) ) CALL iom_put( "dmsspr"      , - wfx_spr                  )                            ! Snow mass change through snow fall 
    280       IF ( iom_use( "dmsssi"   ) ) CALL iom_put( "dmsssi"      ,   wfx_sni*rhosn/rhoic      )                            ! Snow mass change through snow-to-ice conversion 
     280      IF ( iom_use( "dmsssi"   ) ) CALL iom_put( "dmsssi"      ,   wfx_sni*rhosn*r1_rhoic   )                            ! Snow mass change through snow-to-ice conversion 
    281281 
    282282      IF ( iom_use( "dmsmel"   ) ) CALL iom_put( "dmsmel"      , - wfx_snw_sum              )                            ! Snow mass change through melt 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r8233 r8522  
    9191   REAL(wp), PUBLIC ::   r1_rhoic                    !: 1 / rhoic 
    9292   REAL(wp), PUBLIC ::   r1_rhosn                    !: 1 / rhosn 
     93   REAL(wp), PUBLIC ::   r1_cpic                     !: 1 / cpic 
    9394#endif 
    9495   !!---------------------------------------------------------------------- 
     
    159160      r1_rhoic = 1._wp / rhoic 
    160161      r1_rhosn = 1._wp / rhosn 
     162      r1_cpic  = 1._wp / cpic 
    161163#endif 
    162164      IF(lwp) THEN 
Note: See TracChangeset for help on using the changeset viewer.