Changeset 12371


Ignore:
Timestamp:
2020-02-12T13:13:26+01:00 (8 weeks ago)
Author:
dancopsey
Message:

Add changes to match CICE settings.

Location:
NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90

    r12369 r12371  
    127127      REAL(wp), PARAMETER ::   zrmin       = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding 
    128128      REAL(wp), PARAMETER ::   zrmax       = 0.70_wp  ! maximum     -           -         -         -            - 
    129       REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio 
     129      REAL(wp), PARAMETER ::   zpnd_aspect = 0.174_wp   ! pond aspect ratio 
    130130      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature 
    131131      ! 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_zdf_bl99.F90

    r12369 r12371  
    3131   !!---------------------------------------------------------------------- 
    3232   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    33    !! $Id$ 
     33   !! $Id: icethd_zdf_bl99.F90 10926 2019-05-03 12:32:10Z clem $ 
    3434   !! Software governed by the CeCILL license (see ./LICENSE) 
    3535   !!---------------------------------------------------------------------- 
     
    8989      REAL(wp) ::   zg1       =  2._wp        ! 
    9090      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat 
    91       REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13) 
     91      REAL(wp) ::   zbeta     =  0.13_wp     ! for thermal conductivity (could be 0.13) 
    9292      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow 
    9393      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity 
    9494      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C  
    9595      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature  
    96       REAL(wp) ::   zhs_min   =  0.01_wp      ! minimum snow thickness for conductivity calculation  
     96      REAL(wp) ::   zhs_min   =  0.1_wp       ! minimum snow thickness for conductivity calculation  
    9797      REAL(wp) ::   ztmelts                   ! ice melting temperature 
    9898      REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
     
    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)  
     104      LOGICAL,  DIMENSION(jpij) ::   l_snow       ! is snow above the threshold 
    104105      REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration 
    105106      REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness 
     
    153154      !------------------ 
    154155      DO ji = 1, npti 
    155          isnow(ji) = 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) )  ! is there snow or not 
     156          
     157         ! If the snow thickness drops below zhs_min then reduce the snow fraction instead  
     158         IF( h_s_1d(ji) < zhs_min ) THEN  
     159           ! isnow(ji) = h_s_1d(ji) / zhs_min  
     160           isnow(ji) = 0.0_wp 
     161           zh_s(ji) = zhs_min * r1_nlay_s  
     162           l_snow(ji) = .false. 
     163         ELSE  
     164           isnow(ji) = 1.0_wp  
     165           zh_s(ji) = h_s_1d(ji) * r1_nlay_s  
     166           l_snow(ji) = .true. 
     167         END IF  
     168 
    156169         ! layer thickness 
    157170         zh_i(ji) = h_i_1d(ji) * r1_nlay_i 
    158          zh_s(ji) = h_s_1d(ji) * r1_nlay_s 
     171 
    159172         IF( to_print(ji) == 10 ) THEN 
    160173           write(numout,*)'icethd_zdf_bl99: v_i_1d(ji), a_i_1d(ji), h_i_1d(ji) = ',v_i_1d(ji), ' ', a_i_1d(ji), ' ', h_i_1d(ji) 
     
    166179      END WHERE 
    167180      ! 
    168       WHERE( zh_s(1:npti) > 0._wp   )       zh_s(1:npti) = MAX( zhs_min * r1_nlay_s, zh_s(1:npti) ) 
    169181      ! 
    170182      WHERE( zh_s(1:npti) > 0._wp   )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 
     
    310322               IF ( .NOT. l_T_converged(ji) ) & 
    311323                  zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
     324               IF ( .NOT. l_snow(ji) ) zkappa_s(ji,jk) = 0.0_wp 
    312325            END DO 
    313326         END DO 
     
    333346         END DO 
    334347         DO ji = 1, npti   ! Snow-ice interface 
    335             IF ( .NOT. l_T_converged(ji) ) & 
    336                zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
    337348             
    338349            IF (to_print(ji) == 10) THEN 
     
    359370         DO jk = 1, nlay_s 
    360371            DO ji = 1, npti 
    361                zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji) 
     372               IF (l_snow(ji)) then 
     373                  zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji) 
     374               ELSE 
     375                  zeta_s(ji,jk) = 0.0_wp 
     376               END IF 
    362377            END DO 
    363378         END DO 
     
    707722            DO ji = 1, npti 
    708723               !                               !---------------------! 
    709                IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells ! 
     724               IF( l_snow(ji) ) THEN   !  snow-covered cells ! 
    710725                  !                            !---------------------! 
    711726                  ! snow interior terms (bottom equation has the same form as the others) 
     
    895910      ! 
    896911      ! --- calculate conduction fluxes (positive downward) 
    897  
     912      !     bottom ice conduction flux 
    898913      DO ji = 1, npti 
    899          !                                ! surface ice conduction flux 
    900          qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0)      * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  & 
    901             &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0)      * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
    902          !                                ! bottom ice conduction flux 
    903          qcn_ice_bot_1d(ji) =                          - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 
     914         qcn_ice_bot_1d(ji) =  - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 
     915 
     916         IF (to_print(ji) == 1) THEN 
     917           write(numout,*) 'icethd_zdf_bl99: qcn_ice_bot_1d(ji), zkappa_i(ji,nlay_i), zg1, t_bo_1d(ji ), t_i_1d (ji,nlay_i) = ',qcn_ice_bot_1d(ji), zkappa_i(ji,nlay_i), zg1, t_bo_1d(ji ), t_i_1d (ji,nlay_i) 
     918         ENDIF 
    904919      END DO 
    905        
     920      !     surface ice conduction flux 
     921      IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN 
     922         ! 
     923         DO ji = 1, npti 
     924            qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  & 
     925               &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
     926         END DO 
     927         ! 
     928      ELSEIF( k_cnd == np_cnd_ON ) THEN 
     929         ! 
     930         DO ji = 1, npti 
     931            qcn_ice_top_1d(ji) = qcn_ice_1d(ji) 
     932            ! 
     933            t_su_1d(ji) = (  qcn_ice_top_1d(ji) &            ! calculate surface temperature 
     934               &           +           isnow(ji)   * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) & 
     935               &           + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) & 
     936               &          ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 ) 
     937            t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su 
     938         END DO 
     939         ! 
     940      ENDIF 
    906941      ! 
    907942      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! 
     
    911946         DO ji = 1, npti 
    912947            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji)  
    913          END DO 
    914          ! 
    915       ELSEIF( k_cnd == np_cnd_ON ) THEN 
    916          ! 
    917          DO ji = 1, npti 
    918             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qcn_ice_top_1d(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji)  
    919948         END DO 
    920949         ! 
     
    9851014      ! effective conductivity and 1st layer temperature (needed by Met Office) 
    9861015      DO ji = 1, npti 
    987          IF( h_s_1d(ji) > 0.1_wp ) THEN  
    988             cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0) 
     1016         IF( h_i_1d(ji) > 0.1_wp ) THEN 
     1017            cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 
    9891018         ELSE 
    990             IF( h_i_1d(ji) > 0.1_wp ) THEN 
    991                cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 
    992             ELSE 
    993                cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp 
    994             ENDIF 
     1019            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp 
    9951020         ENDIF 
    9961021         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/iceupdate.F90

    r10888 r12371  
    279279      IF( iom_use('hfxcndbot'  ) )   CALL iom_put( "hfxcndbot"  , SUM( qcn_ice_bot * a_i_b, dim=3 ) )   ! Bottom conduction flux 
    280280      IF( iom_use('hfxcndtop'  ) )   CALL iom_put( "hfxcndtop"  , SUM( qcn_ice_top * a_i_b, dim=3 ) )   ! Surface conduction flux 
     281      IF( iom_use('hfxcndcpl'  ) )   CALL iom_put( "hfxcndcpl"  , SUM( qcn_ice * a_i_b, dim=3 ) )       ! Conduction flux we are giving it 
    281282 
    282283      ! diags 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icevar.F90

    r12369 r12371  
    317317      REAL(wp), PARAMETER :: zsi0 = 3.5_wp 
    318318      REAL(wp), PARAMETER :: zsi1 = 4.5_wp 
     319 
     320      REAL(wp), PARAMETER :: min_salin = 0.1_wp     ! Minimum salinity in sea ice 
     321      REAL(wp) ::            zn                     ! Fraction of ice depth 
     322      REAL(wp), PARAMETER :: msal    = 0.573_wp 
     323      REAL(wp), PARAMETER :: nsal    = 0.407_wp 
     324      REAL(wp), PARAMETER :: saltmax = 9.6_wp 
     325      ! REAL(wp), PARAMETER, DIMENSION(4) :: weights = (/0.06363419, 0.25545967, 0.33155539, 0.34935076/) 
     326      REAL(wp), PARAMETER, DIMENSION(4) :: weights = (/0.0703805, 0.25526203, 0.32860314, 0.34575433/) 
     327       
    319328      !!------------------------------------------------------------------- 
    320329 
     
    350359            DO jj = 1, jpj 
    351360               DO ji = 1, jpi 
    352                   zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  ) 
     361                  zalpha(ji,jj,jl) = 1._wp 
    353362                  !                             ! force a constant profile when SSS too low (Baltic Sea) 
    354363                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp   
     
    363372               DO jj = 1, jpj 
    364373                  DO ji = 1, jpi 
    365                      !                          ! linear profile with 0 surface value 
    366                      zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i 
     374 
     375                     !  Copied over from CICE 
     376                     ! zn = (real(jk,kind=dbl_kind)-0.5_wp) / real(nlay_i,kind=dbl_kind) 
     377                     ! sz_i(ji,jj,jk,jl)=(saltmax/2.0_wp)*(1.0_wp-cos(pi*zn**(nsal/(msal+zn)))) 
     378                     ! sz_i(ji,jj,jk,jl) = max(sz_i(ji,jj,jk,jl), min_salin) 
     379 
     380                     ! Weighting method to match CICE 
     381                     zs0 = (s_i(ji,jj,jl)*4.0_wp) * weights(jk) 
    367382                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile 
    368383                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) ) 
     384 
     385 
     386                     ! zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i 
     387                     ! zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile 
     388                     ! sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) ) 
    369389 
    370390                     IF ( jl == 1 .AND. jj == 26 .AND. ji == 42 ) THEN 
    371391                       write(numout,*) 'icevar: jk, sz_i(ji,jj,jk,jl), s_i(ji,jj,jl), weights(jk) = ',jk, ' ',sz_i(ji,jj,jk,jl), ' ',s_i(ji,jj,jl), ' ',weights(jk) 
    372392                     ENDIF 
     393 
    373394                  END DO 
    374395               END DO 
  • NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icewri.F90

    r11081 r12371  
    155155      IF( iom_use('iceconc_cat' ) )   CALL iom_put( "iceconc_cat" , a_i * zmsk00l                                              )   ! area for categories 
    156156      IF( iom_use('icethic_cat' ) )   CALL iom_put( "icethic_cat" , h_i * zmsk00l                                              )   ! thickness for categories 
     157      IF( iom_use('icevol_cat'  ) )   CALL iom_put( "icevol_cat" , v_i * zmsk00l                                               )   ! volume for categories 
    157158      IF( iom_use('snwthic_cat' ) )   CALL iom_put( "snwthic_cat" , h_s * zmsksnl                                              )   ! snow depth for categories 
    158159      IF( iom_use('icesalt_cat' ) )   CALL iom_put( "icesalt_cat" , s_i * zmsk00l                                              )   ! salinity for categories 
Note: See TracChangeset for help on using the changeset viewer.