Changeset 13148


Ignore:
Timestamp:
2020-06-23T22:31:14+02:00 (6 weeks ago)
Author:
clem
Message:

debug heat diffusion routine, so that conduction at the ice surface is correct

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90

    r12919 r13148  
    9595      REAL(wp) ::   zhs_ssl   =  0.03_wp      ! surface scattering layer in the snow  
    9696      REAL(wp) ::   zhi_ssl   =  0.10_wp      ! surface scattering layer in the ice 
     97      REAL(wp) ::   zh_min    =  1.e-3_wp     ! minimum ice/snow thickness for conduction 
    9798      REAL(wp) ::   ztmelts                   ! ice melting temperature 
    9899      REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
     
    123124      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow 
    124125      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow 
     126      REAL(wp), DIMENSION(jpij)            ::   zkappa_comb ! Combined snow and ice surface conductivity 
    125127      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term 
    126128      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term 
     
    130132      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
    131133      REAL(wp), DIMENSION(jpij)            ::   za_s_fra    ! ice fraction covered by snow  
     134      REAL(wp), DIMENSION(jpij)            ::   isnow       ! snow presence (1) or not (0)  
    132135      ! 
    133136      ! Mono-category 
     
    143146      END DO 
    144147 
    145       CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) )   ! calculate ice fraction covered by snow 
     148      ! calculate ice fraction covered by snow for radiation 
     149      CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) ) 
    146150       
    147151      !------------------ 
     
    158162      ENDIF 
    159163      ! 
    160       ! layers thicknesses 
     164      ! thicknesses 
    161165      DO ji = 1, npti 
    162          zh_s(ji) = MAX( zhs_ssl , h_s_1d(ji) ) * r1_nlay_s ! set a minimum snw thickness for conduction 
    163          zh_i(ji) = MAX( zhi_ssl , h_i_1d(ji) ) * r1_nlay_i ! set a minimum ice thickness for conduction 
     166         ! ice thickness 
     167         IF( h_i_1d(ji) > 0._wp ) THEN  
     168            zh_i  (ji) = MAX( zh_min , h_i_1d(ji) ) * r1_nlay_i ! set a minimum thickness for conduction 
     169            z1_h_i(ji) = 1._wp / zh_i(ji)                       !       it must be very small 
     170         ELSE 
     171            zh_i  (ji) = 0._wp 
     172            z1_h_i(ji) = 0._wp 
     173         ENDIF 
     174         ! snow thickness 
     175         IF( h_s_1d(ji) > 0._wp ) THEN 
     176            zh_s  (ji) = MAX( zh_min , h_s_1d(ji) ) * r1_nlay_s ! set a minimum thickness for conduction 
     177            z1_h_s(ji) = 1._wp / zh_s(ji)                       !       it must be very small 
     178            isnow (ji) = 1._wp 
     179         ELSE 
     180            zh_s  (ji) = 0._wp 
     181            z1_h_s(ji) = 0._wp 
     182            isnow (ji) = 0._wp 
     183         ENDIF 
    164184      END DO 
    165       ! 
    166       WHERE( h_i_1d(1:npti) >= zhi_ssl )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti) 
    167       ELSEWHERE                            ;   z1_h_i(1:npti) = 0._wp ! put 0 if ice thick < scattering layer 
    168       END WHERE 
    169       ! 
    170       WHERE( h_s_1d(1:npti) >= zhs_ssl )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 
    171       ELSEWHERE                            ;   z1_h_s(1:npti) = 0._wp ! put 0 if snw thick < scattering layer 
    172       END WHERE 
     185      ! clem: we should apply correction on snow thickness to take into account snow fraction 
     186      !       it must be a distribution, so it is a bit complicated 
    173187      ! 
    174188      ! Store initial temperatures and non solar heat fluxes 
    175189      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN 
    176          ! 
    177190         ztsub      (1:npti) = t_su_1d(1:npti)                          ! surface temperature at iteration n-1 
    178191         ztsuold    (1:npti) = t_su_1d(1:npti)                          ! surface temperature initial value 
     
    200213      END DO 
    201214      ! 
    202       zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * za_s_fra(1:npti) + qtr_ice_top_1d(1:npti) * (1._wp - za_s_fra(1:npti)) 
     215      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * za_s_fra(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - za_s_fra(1:npti) ) 
    203216      DO jk = 1, nlay_i  
    204217         DO ji = 1, npti 
    205218            !                             ! radiation transmitted below the layer-th ice layer 
    206219            zradtr_i(ji,jk) =           za_s_fra(ji)   * zradtr_s(ji,nlay_s)                       &   ! part covered by snow 
    207                &                                       * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) & 
     220               &                                       * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zh_min  ) ) & 
    208221               &            + ( 1._wp - za_s_fra(ji) ) * qtr_ice_top_1d(ji)                        &   ! part snow free 
    209                &                                       * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 
    210              
    211             zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 
     222               &                                       * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) )             
    212223            !                             ! radiation absorbed by the layer-th ice layer 
    213224            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 
     
    304315         END DO 
    305316         DO ji = 1, npti   ! Snow-ice interface 
    306             IF ( .NOT. l_T_converged(ji) ) THEN 
    307                IF( h_s_1d(ji) >= zhs_ssl ) THEN 
    308                   zkappa_s(ji,nlay_s) =   zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / & 
    309                      &                  ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) ) 
    310                ELSE 
    311                   zkappa_s(ji,nlay_s) = 0._wp 
    312                ENDIF 
    313             ENDIF 
     317            IF ( .NOT. l_T_converged(ji) ) & 
     318               zkappa_s(ji,nlay_s) = isnow(ji) * zghe(ji) * rn_cnd_s * ztcond_i(ji,0) & 
     319                  &                            / ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) ) 
    314320         END DO 
    315321 
     
    324330         END DO 
    325331         DO ji = 1, npti   ! Snow-ice interface 
    326             IF ( .NOT. l_T_converged(ji) ) & 
    327                zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * za_s_fra(ji) + zkappa_i(ji,0) * (1._wp - za_s_fra(ji)) 
     332            IF ( .NOT. l_T_converged(ji) ) THEN 
     333               ! Calculate combined surface snow and ice conductivity to pass through the coupler 
     334               zkappa_comb(ji) = isnow(ji) * zkappa_s(ji,0) + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) 
     335               ! If there is snow then use the same snow-ice interface conductivity for the top layer of ice 
     336               IF( h_s_1d(ji) > 0._wp )   zkappa_i(ji,0) = zkappa_s(ji,nlay_s) 
     337           ENDIF 
    328338         END DO 
    329339         ! 
     
    334344            DO ji = 1, npti 
    335345               zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 
    336                zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi )  
     346               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / zcpi 
    337347            END DO 
    338348         END DO 
     
    558568                  ztsub(ji) = t_su_1d(ji) 
    559569                  IF( t_su_1d(ji) < rt0 ) THEN 
    560                      t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    561                         &          ( za_s_fra(ji) * t_s_1d(ji,1) + (1._wp - za_s_fra(ji)) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     570                     t_su_1d(ji) = (  zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
     571                        &           ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
    562572                  ENDIF 
    563573               ENDIF 
    564574            END DO 
     575            !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 
    565576            ! 
    566577            !-------------------------------------------------------------- 
     
    575586 
    576587               IF ( .NOT. l_T_converged(ji) ) THEN 
     588 
    577589                  t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
    578590                  zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 
    579591 
    580                   t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 
    581                   zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 
     592                  IF( h_s_1d(ji) > 0._wp ) THEN 
     593                     DO jk = 1, nlay_s 
     594                        t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
     595                        zdti_max      = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
     596                     END DO 
     597                  ENDIF 
    582598 
    583599                  DO jk = 1, nlay_i 
     
    586602                     zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    587603                  END DO 
     604                   
    588605                  ! convergence test 
    589606                  IF( ln_zdf_chkcvg ) THEN 
     
    745762               ENDIF 
    746763            END DO 
     764            !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 
    747765            ! 
    748766            !-------------------------------------------------------------- 
     
    757775 
    758776               IF ( .NOT. l_T_converged(ji) ) THEN 
    759                   ! t_s 
    760                   t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 
    761                   zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 
    762                   ! t_i 
     777 
     778                  IF( h_s_1d(ji) > 0._wp ) THEN 
     779                     DO jk = 1, nlay_s 
     780                        t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
     781                        zdti_max      = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
     782                     END DO 
     783                  ENDIF 
     784 
    763785                  DO jk = 1, nlay_i 
    764786                     ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
     
    766788                     zdti_max      =  MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    767789                  END DO 
     790 
    768791                  ! convergence test 
    769792                  IF( ln_zdf_chkcvg ) THEN 
     
    789812      !     bottom ice conduction flux 
    790813      DO ji = 1, npti 
    791          qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 
     814         qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 
    792815      END DO 
    793816      !     surface ice conduction flux 
     
    795818         ! 
    796819         DO ji = 1, npti 
    797             qcn_ice_top_1d(ji) =  -           za_s_fra(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) & 
    798                &                  - ( 1._wp - za_s_fra(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
     820            qcn_ice_top_1d(ji) = -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) & 
     821               &                 - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
    799822         END DO 
    800823         ! 
     
    810833         ! 
    811834         DO ji = 1, npti 
    812             t_su_1d(ji) = (  qcn_ice_top_1d(ji) &            ! calculate surface temperature 
    813                &                +          za_s_fra(ji)  * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) & 
    814                &                + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) & 
    815                &          ) / MAX( epsi10, za_s_fra(ji)  * zkappa_s(ji,0) * zg1s + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1 ) 
     835            t_su_1d(ji) = ( qcn_ice_top_1d(ji) +          isnow(ji)   * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) + & 
     836               &                                ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) ) & 
     837               &          / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 ) 
    816838            t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su 
    817839         END DO 
     
    873895      ! this is a conductivity at mid-layer, hence the factor 2 
    874896      DO ji = 1, npti 
    875          IF( h_i_1d(ji) >= zhi_ssl ) THEN   ! zkappa_i=0 when h_i<zhi_ssl (but ztcond_i/=0) 
     897         IF( h_i_1d(ji) >= zhi_ssl ) THEN 
     898!!clem metO            cnd_ice_1d(ji) = 2._wp * zkappa_comb(ji) 
    876899            cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 
    877900         ELSE 
    878901            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / zhi_ssl ! cnd_ice is capped by: cond_i/zhi_ssl 
    879902         ENDIF 
    880          t1_ice_1d(ji) = za_s_fra(ji) * t_s_1d(ji,1) + ( 1._wp - za_s_fra(ji) ) * t_i_1d(ji,1) 
     903         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) 
    881904      END DO 
    882905      ! 
     
    893916         !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    894917         IF( h_s_1d(ji) >= zhs_ssl ) THEN 
    895             t_si_1d(ji) = ( rn_cnd_s * zh_i(ji) * t_s_1d(ji,1) + ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / & 
    896                &          ( rn_cnd_s * zh_i(ji)                + ztcond_i(ji,1) * zh_s(ji)                ) 
     918            t_si_1d(ji) = (   rn_cnd_s       * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,1)   & 
     919               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) & 
     920               &          / ( rn_cnd_s       * h_i_1d(ji) * r1_nlay_i & 
     921               &            + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s ) 
    897922         ELSE 
    898923            t_si_1d(ji) = t_su_1d(ji) 
Note: See TracChangeset for help on using the changeset viewer.