Changeset 5202


Ignore:
Timestamp:
2015-04-07T17:40:16+02:00 (6 years ago)
Author:
clem
Message:

LIM3: important bug fix to avoid crashes. See ticket #1508

Location:
trunk/NEMOGCM/NEMO/LIM_SRC_3
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r5123 r5202  
    314314            DO ji = 1, jpi 
    315315               a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration 
    316                ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))  ! ice thickness 
     316               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))   ! ice thickness 
    317317               ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) )  ! snow depth 
    318                sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin ! salinity 
    319                o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp + ( 1._wp - zswitch(ji,jj) ) ! age 
     318               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity 
     319               o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                        ! age (1 day) 
    320320               t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 
    321321 
     
    333333               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 
    334334               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content 
    335             END DO ! ji 
    336          END DO ! jj 
    337       END DO ! jl 
     335            END DO 
     336         END DO 
     337      END DO 
    338338 
    339339      ! Snow temperature and heat content 
     
    348348                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
    349349                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
    350                END DO ! ji 
    351             END DO ! jj 
    352          END DO ! jl 
    353       END DO ! jk 
     350               END DO 
     351            END DO 
     352         END DO 
     353      END DO 
    354354 
    355355      ! Ice salinity, temperature and heat content 
     
    369369                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
    370370                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
    371                END DO ! ji 
    372             END DO ! jj 
    373          END DO ! jl 
    374       END DO ! jk 
     371               END DO 
     372            END DO 
     373         END DO 
     374      END DO 
    375375 
    376376      tn_ice (:,:,:) = t_su (:,:,:) 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r5181 r5202  
    127127      REAL(wp) ::   za, zfac              ! local scalar 
    128128      CHARACTER (len = 15) ::   fieldid 
    129       REAL(wp), POINTER, DIMENSION(:,:) ::   closing_net     ! net rate at which area is removed    (1/s) 
    130                                                              ! (ridging ice area - area of new ridges) / dt 
    131       REAL(wp), POINTER, DIMENSION(:,:) ::   divu_adv        ! divu as implied by transport scheme  (1/s) 
    132       REAL(wp), POINTER, DIMENSION(:,:) ::   opning          ! rate of opening due to divergence/shear 
    133       REAL(wp), POINTER, DIMENSION(:,:) ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
    134       REAL(wp), POINTER, DIMENSION(:,:) ::   msnow_mlt       ! mass of snow added to ocean (kg m-2) 
    135       REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2) 
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
     129      REAL(wp), POINTER, DIMENSION(:,:)   ::   closing_net     ! net rate at which area is removed    (1/s) 
     130                                                               ! (ridging ice area - area of new ridges) / dt 
     131      REAL(wp), POINTER, DIMENSION(:,:)   ::   divu_adv        ! divu as implied by transport scheme  (1/s) 
     132      REAL(wp), POINTER, DIMENSION(:,:)   ::   opning          ! rate of opening due to divergence/shear 
     133      REAL(wp), POINTER, DIMENSION(:,:)   ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
     134      REAL(wp), POINTER, DIMENSION(:,:)   ::   msnow_mlt       ! mass of snow added to ocean (kg m-2) 
     135      REAL(wp), POINTER, DIMENSION(:,:)   ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2) 
     136      REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
    137137      ! 
    138138      INTEGER, PARAMETER ::   nitermax = 20     
     
    142142      IF( nn_timing == 1 )  CALL timing_start('limitd_me') 
    143143 
    144       CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
     144      CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    145145 
    146146      IF(ln_ctl) THEN 
     
    154154      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    155155 
     156      CALL lim_var_zapsmall 
    156157      CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
     158 
    157159      !-----------------------------------------------------------------------------! 
    158160      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
     
    364366      ENDIF 
    365367 
    366       ! updates 
    367368      CALL lim_var_agg( 1 )  
    368369 
     
    830831      REAL(wp), POINTER, DIMENSION(:,:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges 
    831832      REAL(wp), POINTER, DIMENSION(:,:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice 
    832       REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! areal age content of ridged & rifging ice 
    833833      REAL(wp), POINTER, DIMENSION(:,:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2 
    834834 
     
    839839      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg2   ! sal*volume of new ridges 
    840840      REAL(wp), POINTER, DIMENSION(:,:) ::   smsw    ! sal*volume of water trapped into ridges 
     841      REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! ice age of ice ridged 
    841842 
    842843      REAL(wp), POINTER, DIMENSION(:,:) ::   afrft            ! fraction of category area rafted 
     
    844845      REAL(wp), POINTER, DIMENSION(:,:) ::   virft , vsrft    ! ice & snow volume of rafting ice 
    845846      REAL(wp), POINTER, DIMENSION(:,:) ::   esrft , smrft    ! snow energy & salinity of rafting ice 
    846       REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! areal age content of rafted ice & rafting ice 
     847      REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! ice age of ice rafted 
    847848 
    848849      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eirft      ! ice energy of rafting ice 
     
    854855      CALL wrk_alloc( (jpi+1)*(jpj+1),       indxi, indxj ) 
    855856      CALL wrk_alloc( jpi, jpj,              vice_init, vice_final, eice_init, eice_final ) 
    856       CALL wrk_alloc( jpi, jpj,              afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 ) 
    857       CALL wrk_alloc( jpi, jpj,              vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw ) 
     857      CALL wrk_alloc( jpi, jpj,              afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 
     858      CALL wrk_alloc( jpi, jpj,              vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 
    858859      CALL wrk_alloc( jpi, jpj,              afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    859860      CALL wrk_alloc( jpi, jpj, jpl,         aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
     
    897898      vsnwn_init(:,:,:)   = v_s  (:,:,:) 
    898899      smv_i_init(:,:,:)   = smv_i(:,:,:) 
    899       oa_i_init (:,:,:)   = oa_i (:,:,:) 
    900900      esnwn_init(:,:,:)   = e_s  (:,:,1,:) 
    901901      eicen_init(:,:,:,:) = e_i  (:,:,:,:) 
     902      oa_i_init (:,:,:)   = oa_i (:,:,:) 
    902903 
    903904      ! 
     
    939940            arft2(ji,jj) = arft1(ji,jj) / kraft 
    940941 
    941             oirdg1(ji,jj)= aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice 
    942             oirft1(ji,jj)= araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice 
    943             oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1) 
    944             oirft2(ji,jj)= oirft1(ji,jj) / kraft 
    945  
    946942            !--------------------------------------------------------------- 
    947943            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1  
     
    971967            vsw  (ji,jj) = vrdg1(ji,jj) * rn_por_rdg 
    972968 
    973             vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 
    974             esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 
    975             srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 
     969            vsrdg (ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 
     970            esrdg (ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 
     971            srdg1 (ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 
     972            oirdg1(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) 
     973            oirdg2(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) / krdg(ji,jj,jl1)  
    976974 
    977975            ! rafting volumes, heat contents ... 
    978             virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 
    979             vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 
    980             esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 
    981             smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)  
     976            virft (ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 
     977            vsrft (ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 
     978            esrft (ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 
     979            smrft (ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)  
     980            oirft1(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj)  
     981            oirft2(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) / kraft  
    982982 
    983983            ! substract everything 
    984             a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1(ji,jj)  - arft1(ji,jj) 
    985             v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1(ji,jj)  - virft(ji,jj) 
    986             v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg(ji,jj)  - vsrft(ji,jj) 
    987             e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj)  - esrft(ji,jj) 
     984            a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1 (ji,jj) - arft1 (ji,jj) 
     985            v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1 (ji,jj) - virft (ji,jj) 
     986            v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg (ji,jj) - vsrft (ji,jj) 
     987            e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg (ji,jj) - esrft (ji,jj) 
     988            smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1 (ji,jj) - smrft (ji,jj) 
    988989            oa_i(ji,jj,jl1)  = oa_i(ji,jj,jl1)  - oirdg1(ji,jj) - oirft1(ji,jj) 
    989             smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj)  - smrft(ji,jj) 
    990990 
    991991            !----------------------------------------------------------------- 
    992992            ! 3.5) Compute properties of new ridges 
    993993            !----------------------------------------------------------------- 
    994             !------------- 
     994            !--------- 
    995995            ! Salinity 
    996             !------------- 
     996            !--------- 
    997997            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014 
    998998            srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
     
    11181118               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea 
    11191119 
    1120             END DO ! ij 
     1120            END DO 
    11211121 
    11221122            ! Transfer ice energy to category jl2 by ridging 
     
    11451145                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft 
    11461146                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)     
    1147                   oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)     
     1147                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj) 
    11481148               ENDIF 
    11491149               ! 
     
    11871187      CALL wrk_dealloc( (jpi+1)*(jpj+1),        indxi, indxj ) 
    11881188      CALL wrk_dealloc( jpi, jpj,               vice_init, vice_final, eice_init, eice_final ) 
    1189       CALL wrk_dealloc( jpi, jpj,               afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 ) 
    1190       CALL wrk_dealloc( jpi, jpj,               vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw ) 
     1189      CALL wrk_dealloc( jpi, jpj,               afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 
     1190      CALL wrk_dealloc( jpi, jpj,               vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 
    11911191      CALL wrk_dealloc( jpi, jpj,               afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    11921192      CALL wrk_dealloc( jpi, jpj, jpl,          aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r5134 r5202  
    130130               rswitch           = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) )    !0 if no ice and 1 if yes 
    131131               zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 
    132 !clem               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
    133132               zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl)  
    134133            END DO 
     
    737736      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   ! snow volume summed over categories 
    738737      !!------------------------------------------------------------------ 
    739       !! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate 
    740738       
    741739      CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger 
     
    844842            zdvice(:,:,jl) = 0._wp 
    845843         ENDIF 
    846  
    847 !         ! clem-change begin: why not doing that? 
    848 !         DO jj = 1, jpj 
    849 !            DO ji = 1, jpi 
    850 !               IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 
    851 !                  ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 
    852 !                  a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)  
    853 !               ENDIF 
    854 !            END DO 
    855 !         END DO 
    856          ! clem-change end 
    857844 
    858845      END DO 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5167 r5202  
    8989      REAL(wp) :: zfric_u, zqld, zqfr 
    9090      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    91       REAL(wp), PARAMETER :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04) 
    92       REAL(wp), PARAMETER :: zch        = 0.0057_wp    ! heat transfer coefficient 
     91      REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04) 
     92      REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient 
    9393      ! 
    9494      REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns 
     
    353353      END DO 
    354354  
    355       !------------------------ 
    356       ! Ice natural aging               
    357       !------------------------ 
    358       oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice /rday 
    359  
    360355      !---------------------------------- 
    361356      ! Change thickness to volume 
     
    365360      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
    366361 
     362      ! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting in monocat) 
     363      DO jl  = 1, jpl 
     364         DO jj = 1, jpj 
     365            DO ji = 1, jpi 
     366               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi10 ) ) 
     367               oa_i(ji,jj,jl) = rswitch * oa_i(ji,jj,jl) * a_i(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) 
     368            END DO 
     369         END DO 
     370      END DO 
     371 
    367372      CALL lim_var_zapsmall 
     373 
    368374      !-------------------------------------------- 
    369375      ! Diagnostic thermodynamic growth rates 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5167 r5202  
    9191      REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2) 
    9292      REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2) 
    93       INTEGER , POINTER, DIMENSION(:) ::   icount      ! number of layers vanished by melting  
    9493 
    9594      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt  
     
    9998      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah 
    10099      REAL(wp), POINTER, DIMENSION(:,:) ::   zh_i      ! ice layer thickness 
     100      INTEGER , POINTER, DIMENSION(:,:) ::   icount    ! number of layers vanished by melting  
    101101 
    102102      REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2) 
     
    119119      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
    120120      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
    121       CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i ) 
    122       CALL wrk_alloc( jpij, icount ) 
     121      CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 
     122      CALL wrk_alloc( jpij, nlay_i, icount ) 
    123123       
    124124      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp 
     
    136136      zh_i      (:,:) = 0._wp        
    137137      zdeltah   (:,:) = 0._wp        
    138       icount    (: = 0 
     138      icount    (:,:) = 0 
    139139 
    140140      ! Initialize enthalpy at nlay_i+1 
     
    328328      DO jk = 1, nlay_i 
    329329         DO ji = kideb, kiut 
    330             zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
    331  
    332             ztmelts        = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K] 
    333  
    334             zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
    335  
    336             zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
    337  
    338             IF( zdE < 0._wp ) THEN                                  
    339                zfmdt       = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
    340             ELSE 
    341                zfmdt       = rhoic * zh_i(ji,jk)                   !!! internal melting 
    342                zdE         = 0._wp                                 !   All the layer melts if t_i(jk) = Tmelt (i.e. zdE = 0) 
     330            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K] 
     331             
     332            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
     333 
     334               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]        
     335               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0) 
     336                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted) 
     337               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
     338                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this 
     339               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0 
     340                          
     341               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     342                
     343               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     344 
     345               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
     346               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 
     347                
     348               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     349               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     350                
     351               ! Contribution to mass flux 
     352               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     353 
     354            ELSE                                !!! Surface melting 
     355                
     356               zEi            = - q_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
     357               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
     358               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
     359                
     360               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
     361                
     362               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
     363                
     364               zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     365                
     366               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
     367                
     368               dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     369                
     370               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     371                
     372               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
     373                
     374               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
     375               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     376                
     377               ! Contribution to heat flux [W.m-2], < 0 
     378               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
     379                
     380               ! Total heat flux used in this process [W.m-2], > 0   
     381               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
     382                
     383               ! Contribution to mass flux 
     384               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     385                
    343386            END IF 
    344  
    345             zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
    346  
    347             zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
    348  
    349             zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
    350  
    351             dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
    352  
    353             zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
    354  
    355             zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    356  
    357             ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
    358             sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    359  
    360             ! Contribution to heat flux [W.m-2], < 0 
    361             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    362  
    363             ! Total heat flux used in this process [W.m-2], > 0   
    364             hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    365  
    366             ! Contribution to mass flux 
    367             wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    368             
    369387            ! record which layers have disappeared (for bottom melting)  
    370388            !    => icount=0 : no layer has vanished 
    371389            !    => icount=5 : 5 layers have vanished 
    372             rswitch     = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )  
    373             icount(ji)  = icount(ji) + NINT( rswitch ) 
    374             zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
     390            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )  
     391            icount(ji,jk) = NINT( rswitch ) 
     392            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 
    375393 
    376394            ! update heat content (J.m-2) and layer thickness 
     
    490508      DO jk = nlay_i, 1, -1 
    491509         DO ji = kideb, kiut 
    492             IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji) ) THEN   ! do not calculate where layer has already disappeared by surface melting  
     510            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting  
    493511 
    494512               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K) 
     
    497515 
    498516                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
    499  
    500                   !!zEw               = rcp * ( t_i_1d(ji,jk) - rt0 )  ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 
    501  
    502517                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    503518                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted) 
    504  
    505                   zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing      
    506                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this 
     519                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing      
     520                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this 
    507521 
    508522                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk) 
    509523 
    510                   zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     524                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0 
    511525 
    512526                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
     
    525539               ELSE                               !!! Basal melting 
    526540 
    527                   zEi               = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
    528  
    529                   zEw               = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0) 
    530  
    531                   zdE               = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
    532  
    533                   zfmdt             = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0) 
    534  
    535                   zdeltah(ji,jk)    = - zfmdt * r1_rhoic         ! Gross thickness change 
    536  
    537                   zdeltah(ji,jk)    = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
     541                  zEi             = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
     542                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0) 
     543                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
     544 
     545                  zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0) 
     546 
     547                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change 
     548 
     549                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
    538550                   
    539                   zq_bo(ji)         = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 
    540  
    541                   dh_i_bott(ji)     = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
    542  
    543                   zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
    544  
    545                   zQm               = zfmdt * zEw         ! Heat exchanged with ocean 
     551                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 
     552 
     553                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
     554 
     555                  zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
     556 
     557                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean 
    546558 
    547559                  ! Contribution to heat flux to the ocean [W.m-2], <0   
    548                   hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
     560                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    549561 
    550562                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 
    551                   sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
     563                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 
    552564                   
    553565                  ! Total heat flux used in this process [W.m-2], >0   
    554                   hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
     566                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    555567                   
    556568                  ! Contribution to mass flux 
    557                   wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     569                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    558570 
    559571                  ! update heat content (J.m-2) and layer thickness 
     
    678690      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
    679691      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 
    680       CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) 
    681       CALL wrk_dealloc( jpij, icount ) 
     692      CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 
     693      CALL wrk_dealloc( jpij, nlay_i, icount ) 
    682694      ! 
    683695      ! 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5180 r5202  
    819819      DO jk = 1, nlay_i             ! Sea ice energy of melting 
    820820         DO ji = kideb, kiut 
    821             ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0  
    822             rswitch      = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rt0) - epsi20 ) ) 
    823             q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                                                & 
    824                &                   + lfus * ( 1.0 - rswitch * ( ztmelts-rt0 ) / MIN( t_i_1d(ji,jk) - rt0, -epsi20 ) )   & 
    825                &                   - rcp  *                   ( ztmelts-rt0 )  )  
     821            ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0 
     822            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point 
     823                                                          !   (sometimes dif scheme produces abnormally high temperatures)    
     824            q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                           & 
     825               &                    + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) )   & 
     826               &                    - rcp  *         ( ztmelts-rt0 )  )  
    826827         END DO 
    827828      END DO 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r5182 r5202  
    106106      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d   ! 1-D version of a_i 
    107107      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d   ! 1-D version of v_i 
    108       REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d  ! 1-D version of oa_i 
    109108      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d ! 1-D version of smv_i 
    110109 
     
    119118      CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 
    120119      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 
    121       CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
     120      CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d ) 
    122121      CALL wrk_alloc( jpij,nlay_i,jpl, ze_i_1d ) 
    123122      CALL wrk_alloc( jpi,jpj, zvrel ) 
     
    292291            CALL tab_2d_1d( nbpac, za_i_1d  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    293292            CALL tab_2d_1d( nbpac, zv_i_1d  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    294             CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    295293            CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    296294            DO jk = 1, nlay_i 
    297295               CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 
    298             END DO ! jk 
    299          END DO ! jl 
     296            END DO 
     297         END DO 
    300298 
    301299         CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , jpi, jpj, npac(1:nbpac) ) 
     
    358356         DO ji = 1, nbpac 
    359357            zo_newice(ji) = 0._wp 
    360          END DO ! ji 
     358         END DO 
    361359 
    362360         !------------------- 
     
    480478         ENDDO 
    481479 
    482          !------------ 
    483          ! Update age  
    484          !------------ 
    485          DO jl = 1, jpl 
    486             DO ji = 1, nbpac 
    487                rswitch          = MAX( 0._wp , SIGN( 1._wp , za_i_1d(ji,jl) - epsi20 ) )  ! 0 if no ice and 1 if yes 
    488                zoa_i_1d(ji,jl)  = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch    
    489             END DO  
    490          END DO    
    491  
    492480         !----------------- 
    493481         ! Update salinity 
     
    506494            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 
    507495            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 
    508             CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj ) 
    509496            CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 
    510497            DO jk = 1, nlay_i 
     
    538525      CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 
    539526      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 
    540       CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
     527      CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d ) 
    541528      CALL wrk_dealloc( jpij,nlay_i,jpl, ze_i_1d ) 
    542529      CALL wrk_dealloc( jpi,jpj, zvrel ) 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r5180 r5202  
    150150 
    151151         IF( zcfl > 0.5_wp .AND. lwp )   ncfl = ncfl + 1 
    152          IF( lwp ) THEN 
    153             IF( ncfl > 0 ) THEN    
    154                WRITE(cltmp,'(i6.1)') ncfl 
    155                CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 
    156             ELSE 
    157             !  WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 ' 
    158             ENDIF 
    159          ENDIF 
     152!!         IF( lwp ) THEN 
     153!!            IF( ncfl > 0 ) THEN    
     154!!               WRITE(cltmp,'(i6.1)') ncfl 
     155!!               CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ') 
     156!!            ELSE 
     157!!            !  WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 ' 
     158!!            ENDIF 
     159!!         ENDIF 
    160160 
    161161         !------------------------- 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r5167 r5202  
    8181            DO ji = 1, jpi 
    8282               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    83                   a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     83                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     84                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
    8485               ENDIF 
    8586            END DO 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r5167 r5202  
    7777            ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch 
    7878            IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN 
    79                a_i (ji,jj,1) = a_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin 
     79               a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin 
     80               oa_i(ji,jj,1) = oa_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin 
    8081            ENDIF 
    8182         END DO 
     
    9495            DO ji = 1, jpi 
    9596               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    96                   a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     97                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
     98                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 
    9799               ENDIF 
    98100            END DO 
     
    154156      ! ------------------------------------------------- 
    155157      DO jl  = 1, jpl 
     158         oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice / rday   ! ice natural aging 
    156159         afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 
    157160      END DO 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r5179 r5202  
    175175                  rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes 
    176176                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch 
     177                  !                                      ! bounding salinity 
     178                  sm_i(ji,jj,jl) = MAX( sm_i(ji,jj,jl), rn_simin ) 
    177179               END DO 
    178180            END DO 
     
    199201                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 
    200202                  t_i(ji,jj,jk,jl) = rt0 + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) 
    201                   t_i(ji,jj,jk,jl) = MIN( rt0, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) )       ! 100-rt0 < t_i < rt0 
     203                  t_i(ji,jj,jk,jl) = MIN( ztmelts, MAX( rt0 - 100._wp, t_i(ji,jj,jk,jl) ) )  ! -100 < t_i < ztmelts 
    202204               END DO 
    203205            END DO 
     
    219221                  ! 
    220222                  t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * zq_s + zfac2 ) 
    221                   t_s(ji,jj,jk,jl) = MIN( rt0, MAX( 173.15, t_s(ji,jj,jk,jl) ) )     ! 100-rt0 < t_i < rt0 
     223                  t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) )     ! -100 < t_s < rt0 
    222224               END DO 
    223225            END DO 
     
    264266      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:) 
    265267      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
    266       oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:) 
    267268      ! 
    268269   END SUBROUTINE lim_var_eqv2glo 
     
    346347                     !                                      ! weighting the profile 
    347348                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) 
     349                     !                                      ! bounding salinity 
     350                     s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) ) 
    348351                  END DO 
    349352               END DO 
     
    503506               ! weighting the profile 
    504507               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 
     508               ! bounding salinity 
     509               s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) ) 
    505510            END DO  
    506511         END DO  
     
    556561                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
    557562                  rswitch          = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch 
     563                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 
     564                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  & 
     565                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 
    558566                  zei              = e_i(ji,jj,jk,jl) 
    559567                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch 
     
    569577               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
    570578               rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch 
    571                 
     579               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 
     580               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  & 
     581                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 
    572582               zsal = smv_i(ji,jj,  jl) 
    573583               zvi  = v_i  (ji,jj,  jl) 
  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r5167 r5202  
    7272      ! Mean category values 
    7373      !----------------------------- 
     74      z1_365 = 1._wp / 365._wp 
    7475 
    7576      CALL lim_var_icetm      ! mean sea ice temperature 
     
    112113         CALL lbc_lnk( z2da, 'T', -1. ) 
    113114         CALL lbc_lnk( z2db, 'T', -1. ) 
    114          CALL iom_put( "uice_ipa"     , z2da                )       ! ice velocity u component 
    115          CALL iom_put( "vice_ipa"     , z2db                )       ! ice velocity v component 
     115         CALL iom_put( "uice_ipa"     , z2da             )       ! ice velocity u component 
     116         CALL iom_put( "vice_ipa"     , z2db             )       ! ice velocity v component 
    116117         DO jj = 1, jpj                                  
    117118            DO ji = 1, jpi 
     
    119120            END DO 
    120121         END DO 
    121          CALL iom_put( "icevel"       , z2d                 )       ! ice velocity module 
     122         CALL iom_put( "icevel"       , z2d              )       ! ice velocity module 
    122123      ENDIF 
    123124      ! 
     
    127128            DO jj = 1, jpj 
    128129               DO ji = 1, jpi 
    129                   z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * oa_i(ji,jj,jl) 
     130                  rswitch    = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) ) 
     131                  z2d(ji,jj) = z2d(ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj), 0.1 ) 
    130132               END DO 
    131133            END DO 
    132134         END DO 
    133          z1_365 = 1._wp / 365._wp 
    134          CALL iom_put( "miceage"     , z2d * z1_365         )        ! mean ice age 
     135         CALL iom_put( "miceage"     , z2d * z1_365      )        ! mean ice age 
    135136      ENDIF 
    136137 
     
    141142            END DO 
    142143         END DO 
    143          CALL iom_put( "micet"       , z2d                  )        ! mean ice temperature 
     144         CALL iom_put( "micet"       , z2d               )        ! mean ice temperature 
    144145      ENDIF 
    145146      ! 
     
    153154            END DO 
    154155         END DO 
    155          CALL iom_put( "icest"       , z2d                 )        ! ice surface temperature 
     156         CALL iom_put( "icest"       , z2d              )        ! ice surface temperature 
    156157      ENDIF 
    157158 
     
    163164            END DO 
    164165         END DO 
    165          CALL iom_put( "icecolf"     , z2d                 )        ! frazil ice collection thickness 
     166         CALL iom_put( "icecolf"     , z2d              )        ! frazil ice collection thickness 
    166167      ENDIF 
    167168 
     
    248249            DO jj = 1, jpj 
    249250               DO ji = 1, jpi 
    250                   rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 
    251                   zoi(ji,jj,jl) = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , epsi06 ) * rswitch 
     251                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) ) 
     252                  rswitch = rswitch * MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - 0.1 ) ) 
     253                  zoi(ji,jj,jl) = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , 0.1 ) * rswitch 
    252254               END DO 
    253255            END DO 
    254256         END DO 
    255          CALL iom_put( "iceage_cat"     , zoi        )        ! ice age for categories 
     257         CALL iom_put( "iceage_cat"   , zoi * z1_365 )        ! ice age for categories 
    256258      ENDIF 
    257259 
     
    264266                  DO ji = 1, jpi 
    265267                     rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 
    266                      zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0* & 
     268                     zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0 * & 
    267269                        ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rt0 ), - epsi06 ) ) * & 
    268270                        rswitch * r1_nlay_i 
     
    271273            END DO 
    272274         END DO 
    273          CALL iom_put( "brinevol_cat"     , zei         )        ! brine volume for categories 
     275         CALL iom_put( "brinevol_cat"     , zei      )        ! brine volume for categories 
    274276      ENDIF 
    275277 
Note: See TracChangeset for help on using the changeset viewer.