# Changeset 5079

Ignore:
Timestamp:
2015-02-11T16:32:03+01:00 (6 years ago)
Message:

LIM3: solve ticket #1421

Location:
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO
Files:
13 edited

Unmodified
Removed
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90

 r5070 DO jj = 1, jpj DO ji = 1, jpi ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rt0 IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

 r5070 zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! area zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) )       ! mean salt content zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * e12t(:,:) * tmask(:,:,1) )  ! mean temp content zbg_tem = glob_sum( ( tm_i(:,:) - rt0 ) * vt_i(:,:) * e12t(:,:) * tmask(:,:,1) )  ! mean temp content !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

 r5078 ! surface temperature DO jl = 1, jpl ! loop over categories t_su  (:,:,jl) = rtt * tmask(:,:,1) tn_ice(:,:,jl) = rtt * tmask(:,:,1) t_su  (:,:,jl) = rt0 * tmask(:,:,1) tn_ice(:,:,jl) = rt0 * tmask(:,:,1) END DO sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin ! salinity o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp + ( 1._wp - zswitch(ji,jj) ) ! age t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rtt ! surf temp t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp ! This case below should not be used if (ht_s/ht_i) is ok in namelist DO jj = 1, jpj DO ji = 1, jpi t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rtt t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! Snow energy of melting e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 DO jj = 1, jpj DO ji = 1, jpi t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rtt t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K ! heat content per unit volume e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & +   lfus    * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & -   rcp     * ( ztmelts - rtt ) ) +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & -   rcp     * ( ztmelts - rt0 ) ) ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 DO jl = 1, jpl DO jk = 1, nlay_i t_i(:,:,jk,jl) = rtt * tmask(:,:,1) t_i(:,:,jk,jl) = rt0 * tmask(:,:,1) END DO DO jk = 1, nlay_s t_s(:,:,jk,jl) = rtt * tmask(:,:,1) t_s(:,:,jk,jl) = rt0 * tmask(:,:,1) END DO END DO
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

 r5070 ELSE ht_i(ji,jj,jl)  = 0._wp t_su(ji,jj,jl)  = rtt t_su(ji,jj,jl)  = rt0 ENDIF END DO
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

 r5078 &    (   zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj)               &   ! pfrld already included in coupled mode &    + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)  *     &   ! heat content of precip &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )   & &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )   & &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ) ELSE zqld =  tmask(ji,jj,1) * rdt_ice *  & &      ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) )    & &    + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)  *             &  ! heat content of precip &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )           & &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )           & &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ) ENDIF &    +             pfrld(ji,jj)   * ( zqns(ji,jj) + zqsr(ji,jj) )                                                  & ! latent heat of precip (note that precip is included in qns but not in qns_ice) &    +   ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  & &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) &    +   ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  & &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) ! ----------------------------------------------------------------------------- ! latent heat of precip (note that precip is included in qns but not in qns_ice) &    +      ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)       & &         * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  & &    +      ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt )       & &         * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  & &    +      ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )       & ! heat flux taken from the ocean where there is open water ice formation &    -      qlead(ji,jj) * r1_rdtice                                                                               & DO jk = 1, nlay_i DO ji = kideb, kiut ztmelts       =  -tmut * s_i_1d(ji,jk) + rtt ztmelts       =  -tmut * s_i_1d(ji,jk) + rt0 ! Conversion q(S,T) -> T (second order equation) zaaa          =  cpic zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) * r1_rhoic - lfus zccc          =  lfus * ( ztmelts - rtt ) zbbb          =  ( rcp - cpic ) * ( ztmelts - rt0 ) + q_i_1d(ji,jk) * r1_rhoic - lfus zccc          =  lfus * ( ztmelts - rt0 ) zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) t_i_1d(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) t_i_1d(ji,jk) =  rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) ! mask temperature rswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) t_i_1d(ji,jk) =  rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rtt t_i_1d(ji,jk) =  rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0 END DO END DO
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

 r5078 DO ji = kideb, kiut rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) ztmelts       = rswitch * rtt + ( 1._wp - rswitch ) * rtt ztmelts       = rswitch * rt0 + ( 1._wp - rswitch ) * rt0 zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) !------------------------------------------------------------------------------! DO ji = kideb, kiut IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting ! Contribution to heat flux to the ocean [W.m-2], < 0 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice ht_s_1d(ji)   = 0._wp q_s_1d (ji,1) = 0._wp t_s_1d (ji,1) = rtt t_s_1d (ji,1) = rt0 END IF END DO zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) zqprec   (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) zqprec   (ji) = rhosn * ( cpic * ( rt0 - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp ! heat flux from snow precip (>0, W.m-2) q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) *             & &            ( (   MAX( 0._wp, dh_s_tot(ji) )              ) * zqprec(ji) +  & &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) ) &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk) END DO zEi            = - q_i_1d(ji,jk) * r1_rhoic             ! Specific enthalpy of layer k [J/kg, <0] ztmelts        = - tmut * s_i_1d(ji,jk) + rtt           ! Melting point of layer k [K] ztmelts        = - tmut * s_i_1d(ji,jk) + rt0           ! Melting point of layer k [K] zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] + ( 1. - zswitch_sal ) * sm_i_1d(ji) ! New ice growth ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0) &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   & &               + rcp  * ( ztmelts-rtt ) &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & &               + rcp  * ( ztmelts-rt0 ) zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0) ztmelts        = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0) &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   & &               + rcp  * ( ztmelts-rtt ) &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & &               + rcp  * ( ztmelts-rt0 ) zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) IF(  zf_tt(ji)  >=  0._wp  .AND. jk > icount(ji) ) THEN   ! do not calculate where layer has already disappeared from surface melting ztmelts = - tmut * s_i_1d(ji,jk) + rtt  ! Melting point of layer jk (K) ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K) IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) !!zEw               = rcp * ( t_i_1d(ji,jk) - rtt )  ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) !!zEw               = rcp * ( t_i_1d(ji,jk) - rt0 )  ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) zEi               = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) zEw               = rcp * ( ztmelts - rtt )    ! Specific enthalpy of meltwater (J/kg, <0) zEw               = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0) zdE               = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) DO ji = kideb, kiut rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rtt t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 END DO  ! ji q_s_1d(ji,jk) = ( 1.0 - rswitch ) * q_s_1d(ji,jk) ! recalculate t_s_1d from q_s_1d t_s_1d(ji,jk) = rtt + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) t_s_1d(ji,jk) = rt0 + ( 1._wp - rswitch ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) END DO END DO
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

 r5078 ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr. ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter t_su_1d(ji) =  MIN( t_su_1d(ji), rtt - ztsu_err )       ! necessary t_su_1d(ji) =  MIN( t_su_1d(ji), rt0 - ztsu_err )       ! necessary zerrit (ji) =  1000._wp                                 ! initial value of error END DO IF( nn_ice_thcon == 0 ) THEN      ! Untersteiner (1964) formula DO ji = kideb , kiut ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rtt ) ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) END DO DO ji = kideb , kiut ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  & MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rtt) MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) END DO IF( nn_ice_thcon == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T DO ji = kideb , kiut ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rtt )   & &                   - 0.011_wp * ( t_i_1d(ji,1) - rtt ) ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   & &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) END DO ztcond_i(ji,jk) = rcdic +                                                                       & &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              & &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rtt )   & &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rtt ) &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   & &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) END DO END DO DO ji = kideb , kiut ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rtt )   & &                        - 0.011_wp * ( t_bo_1d(ji) - rtt ) ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   & &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) END DO DO ji = kideb , kiut ztitemp(ji,jk)   = t_i_1d(ji,jk) zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rtt ) * ( ztib(ji,jk) - rtt ), epsi10 ) zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 ) END DO ENDIF IF ( t_su_1d(ji) < rtt ) THEN IF ( t_su_1d(ji) < rt0 ) THEN !------------------------------------------------------------------------------| !------------------------------------------------------------------------------| ! IF ( t_su_1d(ji) < rtt ) THEN IF ( t_su_1d(ji) < rt0 ) THEN ! !------------------------------------------------------------------------------| isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) ) ztsubit(ji) = t_su_1d(ji) IF( t_su_1d(ji) < rtt ) & IF( t_su_1d(ji) < rt0 ) & t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  & &             ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) ! zerrit(ji) is a measure of error, it has to be under terr_dif DO ji = kideb , kiut t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rtt ) , 190._wp  ) t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  ) zerrit(ji)  =  ABS( t_su_1d(ji) - ztsubit(ji) ) END DO DO jk  =  1, nlay_s DO ji = kideb , kiut t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rtt ), 190._wp  ) t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rt0 ), 190._wp  ) zerrit(ji)    = MAX( zerrit(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) ) END DO DO jk  =  1, nlay_i DO ji = kideb , kiut ztmelt_i      = -tmut * s_i_1d(ji,jk) + rtt ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp ) zerrit(ji)    =  MAX( zerrit(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) ) zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  & &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) IF( t_su_1d(ji) < rtt ) THEN  ! case T_su < 0degC IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ELSE                          ! case T_su = 0degC !----------------------------------------- DO ji = kideb, kiut IF( t_su_1d(ji) < rtt ) THEN  ! case T_su < 0degC IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   & &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) DO jk = 1, nlay_i             ! Sea ice energy of melting DO ji = kideb, kiut ztmelts      = - tmut  * s_i_1d(ji,jk) + rtt rswitch      = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi20 ) ) ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0 rswitch      = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rt0) - epsi20 ) ) q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                                                & &                   + lfus * ( 1.0 - rswitch * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk) - rtt, -epsi20 ) )   & &                   - rcp  *                   ( ztmelts-rtt )  ) &                   + lfus * ( 1.0 - rswitch * ( ztmelts-rt0 ) / MIN( t_i_1d(ji,jk) - rt0, -epsi20 ) )   & &                   - rcp  *                   ( ztmelts-rt0 )  ) END DO END DO DO jk = 1, nlay_s             ! Snow energy of melting DO ji = kideb, kiut q_s_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) q_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) END DO END DO
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

 r5078 ! We assume that new ice is formed at the seawater freezing point DO ji = 1, nbpac ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K) ztmelts       = - tmut * zs_newice(ji) + rt0                  ! Melting point (K) ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                             & &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) )   & &                       - rcp  *         ( ztmelts - rtt )  ) END DO ! ji &                       + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & &                       - rcp  *         ( ztmelts - rt0 )  ) END DO !----------------
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90

 r5070 ! Switches !---------- iflush  = MAX( 0._wp , SIGN( 1._wp , t_su_1d(ji) - rtt )        )     ! =1 if summer iflush  = MAX( 0._wp , SIGN( 1._wp , t_su_1d(ji) - rt0 )        )     ! =1 if summer igravdr = MAX( 0._wp , SIGN( 1._wp , t_bo_1d(ji) - t_su_1d(ji) ) )    ! =1 if t_su < t_bo
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

 r5078 rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes zq_i    = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt              ! Ice layer melt temperature ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0              ! Ice layer melt temperature ! zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation) zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i * r1_rhoic - lfus zccc       =  lfus * (ztmelts-rtt) zbbb       =  ( rcp - cpic ) * ( ztmelts - rt0 ) + zq_i * r1_rhoic - lfus zccc       =  lfus * (ztmelts-rt0) zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) t_i(ji,jj,jk,jl) = rtt + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) t_i(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) )       ! 100-rtt < t_i < rtt t_i(ji,jj,jk,jl) = rt0 + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) t_i(ji,jj,jk,jl) = MIN( rt0, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) )       ! 100-rt0 < t_i < rt0 END DO END DO zq_s  = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp) ! t_s(ji,jj,jk,jl) = rtt + rswitch * ( - zfac1 * zq_s + zfac2 ) t_s(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15, t_s(ji,jj,jk,jl) ) )     ! 100-rtt < t_i < rtt t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * zq_s + zfac2 ) t_s(ji,jj,jk,jl) = MIN( rt0, MAX( 173.15, t_s(ji,jj,jk,jl) ) )     ! 100-rt0 < t_i < rt0 END DO END DO DO jj = 1, jpj DO ji = 1, jpi rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) )  ) zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   & rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  ) zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )   & &                   * v_i(ji,jj,jl) * r1_nlay_i rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) zei              = e_i(ji,jj,jk,jl) e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rtt * ( 1._wp - rswitch ) t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch ) ! update exchanges with ocean hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 ! Zap snow energy !----------------------------------------------------------------- t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rtt * ( 1._wp - rswitch ) t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch ) e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

 r5078 DO jj = 1, jpj DO ji = 1, jpi z2d(ji,jj) = ( tm_i(ji,jj) - rtt ) * zswi(ji,jj) z2d(ji,jj) = ( tm_i(ji,jj) - rt0 ) * zswi(ji,jj) END DO END DO DO jj = 1, jpj DO ji = 1, jpi z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * ( t_su(ji,jj,jl) - rtt ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 ) z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 ) END DO END DO rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0* & ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rtt ), - epsi06 ) ) * & ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rt0 ), - epsi06 ) ) * & rswitch * r1_nlay_i END DO CALL histwrite( kid, "iicethic", kt, icethi        , jpi*jpj, (/1/) ) CALL histwrite( kid, "iiceconc", kt, at_i          , jpi*jpj, (/1/) ) CALL histwrite( kid, "iicetemp", kt, tm_i - rtt    , jpi*jpj, (/1/) ) CALL histwrite( kid, "iicetemp", kt, tm_i - rt0    , jpi*jpj, (/1/) ) CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) ) CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) )
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90

 r5064 t_su(ji,jj,jl)   = rswitch * rn_ice_tem(ib_bdy)  + ( 1.0 - rswitch ) * rn_ice_tem(ib_bdy) DO jk = 1, nlay_s t_s(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rtt t_s(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0 END DO DO jk = 1, nlay_i t_i(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rtt t_i(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0 s_i(ji,jj,jk,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * s_i_min END DO sm_i(ji,jj,jl)   = rswitch * sm_i(ii,ij,jl)  + ( 1.0 - rswitch ) * s_i_min o_i(ji,jj,jl)    = rswitch * o_i(ii,ij,jl)   + ( 1.0 - rswitch ) t_su(ji,jj,jl)   = rswitch * t_su(ii,ij,jl)  + ( 1.0 - rswitch ) * rtt t_su(ji,jj,jl)   = rswitch * t_su(ii,ij,jl)  + ( 1.0 - rswitch ) * rt0 DO jk = 1, nlay_s t_s(ji,jj,jk,jl) = rswitch * t_s(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rtt t_s(ji,jj,jk,jl) = rswitch * t_s(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0 END DO DO jk = 1, nlay_i t_i(ji,jj,jk,jl) = rswitch * t_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rtt t_i(ji,jj,jk,jl) = rswitch * t_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0 s_i(ji,jj,jk,jl) = rswitch * s_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * s_i_min END DO DO jk = 1, nlay_s ! Snow energy of melting e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) ! Multiply by volume, so that heat content in J/m2 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) / nlay_s END DO DO jk = 1, nlay_i ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K ! heat content per unit volume e_i(ji,jj,jk,jl) = rswitch * rhoic * & (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & - rcp      * ( ztmelts - rtt ) ) +   lfus    * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & - rcp      * ( ztmelts - rt0 ) ) ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i
• ## branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

 r5078 REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin] #if defined key_lim3 REAL(wp), PUBLIC ::   rt0_snow = 273.16_wp        !: melting point of snow         [Kelvin] REAL(wp), PUBLIC ::   rt0_ice  = 273.16_wp        !: melting point of ice          [Kelvin] REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin] REAL(wp), PUBLIC ::   rt0_ice  = 273.15_wp        !: melting point of ice          [Kelvin] #else REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin]
Note: See TracChangeset for help on using the changeset viewer.