Changeset 12938


Ignore:
Timestamp:
2020-05-15T18:26:45+02:00 (10 days ago)
Author:
dancopsey
Message:

Merge in missing icethd_zdf_bl99 changes including:

  • Make isnow a float
  • Turn off 0.1m snow test
  • Use conductive heat flux from coupler
File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl_v2/src/ICE/icethd_zdf_bl99.F90

    r11715 r12938  
    101101      REAL(wp) ::   zfac                      ! dummy factor 
    102102      ! 
    103       REAL(wp), DIMENSION(jpij) ::   isnow        ! switch for presence (1) or absence (0) of snow 
     103      REAL(wp), DIMENSION(jpij) ::   isnow        ! fraction of sea ice that is snow covered (for thermodynamic use only) 
    104104      REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration 
    105105      REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness 
     
    147147      !------------------ 
    148148      DO ji = 1, npti 
    149          isnow(ji) = 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) )  ! is there snow or not 
     149 
     150         ! If the snow thickness drops below zhs_min then reduce the snow fraction instead 
     151         IF( h_s_1d(ji) < zhs_min ) THEN 
     152           isnow(ji) = h_s_1d(ji) / zhs_min 
     153           zh_s(ji) = zhs_min * r1_nlay_s 
     154         ELSE 
     155           isnow(ji) = 1.0_wp 
     156           zh_s(ji) = h_s_1d(ji) * r1_nlay_s 
     157         END IF 
     158 
    150159         ! layer thickness 
    151160         zh_i(ji) = h_i_1d(ji) * r1_nlay_i 
    152          zh_s(ji) = h_s_1d(ji) * r1_nlay_s 
     161 
    153162      END DO 
    154163      ! 
     
    156165      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp 
    157166      END WHERE 
    158       ! 
    159       WHERE( zh_s(1:npti) > 0._wp   )       zh_s(1:npti) = MAX( zhs_min * r1_nlay_s, zh_s(1:npti) ) 
    160167      ! 
    161168      WHERE( zh_s(1:npti) > 0._wp   )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 
     
    769776      ! 
    770777      ! --- calculate conduction fluxes (positive downward) 
    771  
     778      !     bottom ice conduction flux 
    772779      DO ji = 1, npti 
    773          !                                ! surface ice conduction flux 
    774          qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0)      * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  & 
    775             &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0)      * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
    776          !                                ! bottom ice conduction flux 
    777          qcn_ice_bot_1d(ji) =                          - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 
     780         qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 
    778781      END DO 
    779        
     782      !     surface ice conduction flux 
     783      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN 
     784         ! 
     785         DO ji = 1, npti 
     786            qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  & 
     787               &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
     788         END DO 
     789         ! 
     790      ELSEIF( k_cnd == np_cnd_ON ) THEN 
     791         ! 
     792         DO ji = 1, npti 
     793            qcn_ice_top_1d(ji) = qcn_ice_1d(ji) 
     794         END DO 
     795         ! 
     796      ENDIF 
     797      !     surface ice temperature 
     798      IF( k_cnd == np_cnd_ON .AND. ln_cndemulate ) THEN 
     799         ! 
     800         DO ji = 1, npti 
     801            t_su_1d(ji) = (  qcn_ice_top_1d(ji) &            ! calculate surface temperature 
     802               &           +           isnow(ji)   * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) & 
     803               &           + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) & 
     804               &          ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 ) 
     805            t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su 
     806         END DO 
     807         ! 
     808      ENDIF 
    780809      ! 
    781810      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! 
     
    787816         END DO 
    788817         ! 
    789       ELSEIF( k_cnd == np_cnd_ON ) THEN 
    790          ! 
    791          DO ji = 1, npti 
    792             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qcn_ice_top_1d(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji)  
    793          END DO 
    794          ! 
    795818      ENDIF 
    796        
    797819      ! 
    798820      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- ! 
     
    839861      ! effective conductivity and 1st layer temperature (needed by Met Office) 
    840862      DO ji = 1, npti 
    841          IF( h_s_1d(ji) > 0.1_wp ) THEN  
    842             cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0) 
     863         IF( h_i_1d(ji) > 0.1_wp ) THEN 
     864            cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 
    843865         ELSE 
    844             IF( h_i_1d(ji) > 0.1_wp ) THEN 
    845                cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 
    846             ELSE 
    847                cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp 
    848             ENDIF 
     866            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp 
    849867         ENDIF 
    850868         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) 
     
    856874         t_i_1d    (1:npti,:) = ztiold        (1:npti,:) 
    857875         qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti) 
    858  
    859          !!clem 
    860          ! remettre t_su_1d, qns_ice_1d et dqns_ice_1d comme avant puisqu'on devrait faire comme si on avant conduction = input 
    861          !clem 
    862876      ENDIF 
    863877      ! 
Note: See TracChangeset for help on using the changeset viewer.