Changeset 12890


Ignore:
Timestamp:
2020-05-07T16:11:23+02:00 (7 months ago)
Author:
clem
Message:

1) implement the optional nn_snwfra, i.e. the fraction of ice covered by sea ice. 2) change some namelists parameters names to avoid confusion. 3) correct a bug at ice initialization

Location:
NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/cfgs/SHARED/namelist_ice_ref

    r12832 r12890  
    5656   rn_ishlat        =   2.            !  lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 
    5757   ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016 
    58       rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast 
     58      rn_lf_depfra  =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast 
    5959                                      !          recommended range: [0.1 ; 0.25] 
    60       rn_icebfr     =  15.            !        maximum bottom stress per unit volume [N/m3] 
    61       rn_lfrelax    =   1.e-5         !        relaxation time scale to reach static friction [s-1] 
    62       rn_tensile    =   0.05          !        isotropic tensile strength [0-0.5??] 
     60      rn_lf_bfr     =  15.            !        maximum bottom stress per unit volume [N/m3] 
     61      rn_lf_relax   =   1.e-5         !        relaxation time scale to reach static friction [s-1] 
     62      rn_lf_tensile =   0.05          !        isotropic tensile strength [0-0.5??] 
    6363/ 
    6464!------------------------------------------------------------------------------ 
     
    110110!------------------------------------------------------------------------------ 
    111111   rn_cio           =   5.0e-03       !  ice-ocean drag coefficient (-) 
    112    rn_blow_s        =   0.66          !  mesure of snow blowing into the leads 
     112   nn_snwfra        =   0             !  calculate the fraction of ice covered by snow (for zdf and albedo) 
     113                                      !     = 0  fraction = 1 (if snow) or 0 (if no snow) 
     114                                      !     = 1  fraction = 1-exp(-0.2*rhos*hsnw) [MetO formulation] 
     115                                      !     = 2  fraction = hsnw / (hsnw+0.02)    [CICE formulation] 
     116   rn_snwblow       =   0.66          !  mesure of snow blowing into the leads 
    113117                                      !     = 1 => no snow blowing, < 1 => some snow blowing 
    114118   nn_flxdist       =  -1             !  Redistribute heat flux over ice categories 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/doc/namelists/namdyn

    r11703 r12890  
    1010   rn_ishlat        =   2.            !  lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 
    1111   ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016 
    12       rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast 
     12      rn_lf_depfra  =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast 
    1313                                      !          recommended range: [0.1 ; 0.25] 
    14       rn_icebfr     =  15.            !        maximum bottom stress per unit volume [N/m3] 
    15       rn_lfrelax    =   1.e-5         !        relaxation time scale to reach static friction [s-1] 
    16       rn_tensile    =   0.05          !        isotropic tensile strength [0-0.5??] 
     14      rn_lf_bfr     =  15.            !        maximum bottom stress per unit volume [N/m3] 
     15      rn_lf_relax   =   1.e-5         !        relaxation time scale to reach static friction [s-1] 
     16      rn_lf_tensile =   0.05          !        isotropic tensile strength [0-0.5??] 
    1717/ 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/ice.F90

    r12832 r12890  
    141141   REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice 
    142142   LOGICAL , PUBLIC ::   ln_landfast_L16  !: landfast ice parameterizationfrom lemieux2016  
    143    REAL(wp), PUBLIC ::   rn_depfra        !:    fraction of ocean depth that ice must reach to initiate landfast ice 
    144    REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home)  
    145    REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction 
    146    REAL(wp), PUBLIC ::   rn_tensile       !:    isotropic tensile strength 
     143   REAL(wp), PUBLIC ::   rn_lf_depfra     !:    fraction of ocean depth that ice must reach to initiate landfast ice 
     144   REAL(wp), PUBLIC ::   rn_lf_bfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home)  
     145   REAL(wp), PUBLIC ::   rn_lf_relax      !:    relaxation time scale (s-1) to reach static friction 
     146   REAL(wp), PUBLIC ::   rn_lf_tensile    !:    isotropic tensile strength 
    147147   ! 
    148148   !                                     !!** ice-ridging/rafting namelist (namdyn_rdgrft) ** 
     
    163163   !                                     !!** ice-surface boundary conditions namelist (namsbc) ** 
    164164                                          ! -- icethd_dh -- ! 
    165    REAL(wp), PUBLIC ::   rn_blow_s        !: coef. for partitioning of snowfall between leads and sea ice 
     165   REAL(wp), PUBLIC ::   rn_snwblow       !: coef. for partitioning of snowfall between leads and sea ice 
     166                                          ! -- icethd_zdf and icealb -- ! 
     167   INTEGER , PUBLIC ::   nn_snwfra        !: calculate the fraction of ice covered by snow 
     168   !                                      !   = 0  fraction = 1 (if snow) or 0 (if no snow) 
     169   !                                      !   = 1  fraction = 1-exp(-0.2*rhos*hsnw) [MetO formulation] 
     170   !                                      !   = 2  fraction = hsnw / (hsnw+0.02)    [CICE formulation] 
    166171                                          ! -- icethd -- ! 
    167172   REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn.F90

    r12720 r12890  
    223223      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  & 
    224224         &             rn_ishlat ,                                                           & 
    225          &             ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 
     225         &             ln_landfast_L16, rn_lf_depfra, rn_lf_bfr, rn_lf_relax, rn_lf_tensile 
    226226      !!------------------------------------------------------------------- 
    227227      ! 
     
    246246         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat 
    247247         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16 
    248          WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra 
    249          WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr 
    250          WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax 
    251          WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile 
     248         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_lf_depfra    = ', rn_lf_depfra 
     249         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_lf_bfr       = ', rn_lf_bfr 
     250         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lf_relax     = ', rn_lf_relax 
     251         WRITE(numout,*) '         isotropic tensile strength                          rn_lf_tensile   = ', rn_lf_tensile 
    252252         WRITE(numout,*) 
    253253      ENDIF 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icedyn_rhg_evp.F90

    r12741 r12890  
    258258 
    259259      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
    260       IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile 
     260      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile 
    261261      ELSE                         ;   zkt = 0._wp 
    262262      ENDIF 
     
    329329               zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    330330               ! ice-bottom stress at U points 
    331                zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 
    332                ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
     331               zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj) 
     332               ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    333333               ! ice-bottom stress at V points 
    334                zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 
    335                ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
     334               zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj) 
     335               ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    336336               ! ice_bottom stress at T points 
    337                zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 
    338                tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
     337               zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj) 
     338               tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    339339            END DO 
    340340         END DO 
     
    506506                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    507507                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    508                         &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     508                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )          & ! static friction => slow decrease to v=0 
    509509                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    510510                        &           )   * zmsk00y(ji,jj) 
     
    513513                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    514514                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    515                         &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     515                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
    516516                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    517517                        &            )   * zmsk00y(ji,jj) 
     
    557557                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    558558                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    559                         &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     559                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )          & ! static friction => slow decrease to v=0 
    560560                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    561561                        &           )   * zmsk00x(ji,jj) 
     
    564564                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    565565                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    566                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     566                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
    567567                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    568568                        &            )   * zmsk00x(ji,jj) 
     
    610610                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    611611                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    612                         &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     612                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )          & ! static friction => slow decrease to v=0 
    613613                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    614614                        &           )   * zmsk00x(ji,jj) 
     
    617617                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    618618                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    619                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     619                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
    620620                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    621621                        &            )   * zmsk00x(ji,jj) 
     
    661661                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    662662                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    663                         &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     663                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )          & ! static friction => slow decrease to v=0 
    664664                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    665665                        &           )   * zmsk00y(ji,jj) 
     
    668668                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    669669                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    670                         &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     670                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
    671671                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    672672                        &            )   * zmsk00y(ji,jj) 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/iceistate.F90

    r12816 r12890  
    199199               si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    200200               si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    201             ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2 
    202                si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 ) 
    203             ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2 
    204                si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 ) 
    205             ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_su, set T_su = T_s 
    206                si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1) 
    207             ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_su, set T_su = T_i 
    208                si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
    209             ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_s, set T_s = T_su 
    210                si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1) 
    211             ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_s, set T_s = T_i 
    212                si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
    213201            ENDIF 
     202            IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2 
     203               &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 ) 
     204            IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2 
     205               &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 ) 
     206            IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_su, set T_su = T_s 
     207               &     si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1) 
     208            IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_su, set T_su = T_i 
     209               &     si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
     210            IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_s, set T_s = T_su 
     211               &     si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1) 
     212            IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_s, set T_s = T_i 
     213               &     si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
    214214            ! 
    215215            ! pond concentration 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icesbc.F90

    r12785 r12890  
    276276      INTEGER ::   ios, ioptio   ! Local integer 
    277277      !! 
    278       NAMELIST/namsbc/ rn_cio, rn_blow_s, nn_flxdist, ln_cndflx, ln_cndemulate 
     278      NAMELIST/namsbc/ rn_cio, rn_snwblow, nn_snwfra, nn_flxdist, ln_cndflx, ln_cndemulate 
    279279      !!------------------------------------------------------------------- 
    280280      ! 
     
    293293         WRITE(numout,*) '   Namelist namsbc:' 
    294294         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio        = ', rn_cio 
    295          WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s     = ', rn_blow_s 
     295         WRITE(numout,*) '      fraction of ice covered by snow (options 0,1,2)  nn_snwfra     = ', nn_snwfra 
     296         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_snwblow    = ', rn_snwblow 
    296297         WRITE(numout,*) '      Multicategory heat flux formulation              nn_flxdist    = ', nn_flxdist 
    297298         WRITE(numout,*) '      Use conduction flux as surface condition         ln_cndflx     = ', ln_cndflx 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_dh.F90

    r10786 r12890  
    1313   !!---------------------------------------------------------------------- 
    1414   !!   ice_thd_dh        : vertical sea-ice growth and melt 
    15    !!   ice_thd_snwblow   : distribute snow fall between ice and ocean 
    16   !!---------------------------------------------------------------------- 
     15   !!---------------------------------------------------------------------- 
    1716   USE dom_oce        ! ocean space and time domain 
    1817   USE phycst         ! physical constants 
     
    2019   USE ice1D          ! sea-ice: thermodynamics variables 
    2120   USE icethd_sal     ! sea-ice: salinity profiles 
     21   USE icevar         ! for CALL ice_var_snwblow 
    2222   ! 
    2323   USE in_out_manager ! I/O manager 
     
    2929 
    3030   PUBLIC   ice_thd_dh        ! called by ice_thd 
    31    PUBLIC   ice_thd_snwblow   ! called in sbcblk/sbccpl and here 
    32  
    33    INTERFACE ice_thd_snwblow 
    34       MODULE PROCEDURE ice_thd_snwblow_1d, ice_thd_snwblow_2d 
    35    END INTERFACE 
    3631 
    3732   !!---------------------------------------------------------------------- 
     
    186181      ! Snow precipitation 
    187182      !------------------- 
    188       CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing 
     183      CALL ice_var_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing 
    189184 
    190185      zdeltah(1:npti,:) = 0._wp 
     
    636631   END SUBROUTINE ice_thd_dh 
    637632 
    638  
    639    !!-------------------------------------------------------------------------- 
    640    !! INTERFACE ice_thd_snwblow 
    641    !! 
    642    !! ** Purpose :   Compute distribution of precip over the ice 
    643    !! 
    644    !!                Snow accumulation in one thermodynamic time step 
    645    !!                snowfall is partitionned between leads and ice. 
    646    !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads 
    647    !!                but because of the winds, more snow falls on leads than on sea ice 
    648    !!                and a greater fraction (1-at_i)^beta of the total mass of snow  
    649    !!                (beta < 1) falls in leads. 
    650    !!                In reality, beta depends on wind speed,  
    651    !!                and should decrease with increasing wind speed but here, it is  
    652    !!                considered as a constant. an average value is 0.66 
    653    !!-------------------------------------------------------------------------- 
    654 !!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE.... 
    655    SUBROUTINE ice_thd_snwblow_2d( pin, pout ) 
    656       REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b ) 
    657       REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout 
    658       pout = ( 1._wp - ( pin )**rn_blow_s ) 
    659    END SUBROUTINE ice_thd_snwblow_2d 
    660  
    661    SUBROUTINE ice_thd_snwblow_1d( pin, pout ) 
    662       REAL(wp), DIMENSION(:), INTENT(in   ) :: pin 
    663       REAL(wp), DIMENSION(:), INTENT(inout) :: pout 
    664       pout = ( 1._wp - ( pin )**rn_blow_s ) 
    665    END SUBROUTINE ice_thd_snwblow_1d 
    666  
    667633#else 
    668634   !!---------------------------------------------------------------------- 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_zdf_bl99.F90

    r12854 r12890  
    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.1_wp       ! minimum snow thickness for conductivity calculation  
     96      REAL(wp) ::   zhs_ssl   =  0.03_wp      ! surface scattering layer in the snow  
     97      REAL(wp) ::   zhi_ssl   =  0.10_wp      ! surface scattering layer in the ice 
    9798      REAL(wp) ::   ztmelts                   ! ice melting temperature 
    9899      REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
    99100      REAL(wp) ::   zcpi                      ! Ice specific heat 
    100101      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat 
    101       REAL(wp) ::   zfac                      ! dummy factor 
    102       ! 
    103       REAL(wp), DIMENSION(jpij) ::   isnow        ! fraction of sea ice that is snow covered (for thermodynamic use only) 
     102      ! 
    104103      REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration 
    105104      REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness 
     
    150149      !------------------ 
    151150      DO ji = 1, npti 
    152  
    153          ! If the snow thickness drops below zhs_min then reduce the snow fraction instead 
    154          !!IF( h_s_1d(ji) < zhs_min ) THEN 
    155          !!  isnow(ji) = h_s_1d(ji) / zhs_min 
    156          !!  zh_s(ji) = zhs_min * r1_nlay_s 
    157          !!ELSE 
    158          !!  isnow(ji) = 1.0_wp 
    159          !!  zh_s(ji) = h_s_1d(ji) * r1_nlay_s 
    160          !!END IF 
    161          isnow(ji) = za_s_fra(ji) !!clem: 2 variables for the same thing 
    162          zh_s(ji) = h_s_1d(ji) * r1_nlay_s 
    163           
    164           
    165          ! layer thickness 
    166          zh_i(ji) = h_i_1d(ji) * r1_nlay_i 
    167  
     151         zh_s(ji) = MAX( zhs_ssl , h_s_1d(ji) ) * r1_nlay_s ! set a minimum snw thickness for conduction 
     152         zh_i(ji) = MAX( zhi_ssl , h_i_1d(ji) ) * r1_nlay_i ! set a minimum ice thickness for conduction 
    168153      END DO 
    169154      ! 
    170       WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti) 
    171       ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp 
     155      WHERE( h_i_1d(1:npti) >= zhi_ssl )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti) 
     156      ELSEWHERE                            ;   z1_h_i(1:npti) = 0._wp ! put 0 if ice thick < scattering layer 
    172157      END WHERE 
    173158      ! 
    174       WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 
    175       ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp 
     159      WHERE( h_s_1d(1:npti) >= zhs_ssl )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 
     160      ELSEWHERE                            ;   z1_h_s(1:npti) = 0._wp ! put 0 if snw thick < scattering layer 
    176161      END WHERE 
    177162      ! 
     
    198183         DO ji = 1, npti 
    199184            !                             ! radiation transmitted below the layer-th snow layer 
    200             zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) ) 
     185            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * MAX( 0._wp, zh_s(ji) * REAL(jk) - zhs_ssl ) ) 
    201186            !                             ! radiation absorbed by the layer-th snow layer 
    202187            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 
     
    204189      END DO 
    205190      ! 
    206       zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - isnow(1:npti) ) 
     191      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)) 
    207192      DO jk = 1, nlay_i  
    208193         DO ji = 1, npti 
    209194            !                             ! radiation transmitted below the layer-th ice layer 
    210             zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) 
     195            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) ) 
    211196            !                             ! radiation absorbed by the layer-th ice layer 
    212197            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 
     
    216201      qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice 
    217202      ! 
    218       iconv    = 0          ! number of iterations 
     203      iconv = 0          ! number of iterations 
    219204      ! 
    220205      l_T_converged(:) = .FALSE. 
     
    243228               DO ji = 1, npti 
    244229                  ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  & 
    245                      &                         MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
     230                     &                    MIN( -epsi10, 0.5_wp * (  t_i_1d(ji,jk) +  t_i_1d(ji,jk+1) ) - rt0 ) 
    246231               END DO 
    247232            END DO 
     
    251236            DO ji = 1, npti 
    252237               ztcond_i_cp(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
    253                   &                           - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
     238                  &                            - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
    254239               ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
    255                   &                           - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
     240                  &                            - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
    256241            END DO 
    257242            DO jk = 1, nlay_i-1 
    258243               DO ji = 1, npti 
    259                   ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
    260                      &                        MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) & 
    261                      &                       - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
     244                  ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /       & 
     245                     &                         MIN( -epsi10, 0.5_wp * (  t_i_1d(ji,jk) +  t_i_1d(ji,jk+1) ) - rt0 ) & 
     246                     &                        - 0.011_wp * ( 0.5_wp * (  t_i_1d(ji,jk) +  t_i_1d(ji,jk+1) ) - rt0 ) 
    262247               END DO 
    263248            END DO 
     
    304289         DO ji = 1, npti   ! Snow-ice interface 
    305290            IF ( .NOT. l_T_converged(ji) ) THEN 
    306                zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 
    307                IF( zfac > epsi10 ) THEN 
    308                   zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 
     291               IF( h_s_1d(ji) >= zhs_ssl ) THEN 
     292                  zkappa_s(ji,nlay_s) =   zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / & 
     293                     &                  ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) ) 
    309294               ELSE 
    310295                  zkappa_s(ji,nlay_s) = 0._wp 
     
    324309         DO ji = 1, npti   ! Snow-ice interface 
    325310            IF ( .NOT. l_T_converged(ji) ) & 
    326                zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
     311               zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * za_s_fra(ji) + zkappa_i(ji,0) * (1._wp - za_s_fra(ji)) 
    327312         END DO 
    328313         ! 
     
    558543                  IF( t_su_1d(ji) < rt0 ) THEN 
    559544                     t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    560                         &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     545                        &          ( 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)) 
    561546                  ENDIF 
    562547               ENDIF 
     
    790775         ! 
    791776         DO ji = 1, npti 
    792             qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  & 
    793                &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
     777            qcn_ice_top_1d(ji) =  -           za_s_fra(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  & 
     778               &                  - ( 1._wp - za_s_fra(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
    794779         END DO 
    795780         ! 
     
    806791         DO ji = 1, npti 
    807792            t_su_1d(ji) = (  qcn_ice_top_1d(ji) &            ! calculate surface temperature 
    808                &           +           isnow(ji)   * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) & 
    809                &           + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) & 
    810                &          ) / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 ) 
     793               &                +          za_s_fra(ji)  * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) & 
     794               &                + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1  * t_i_1d(ji,1) & 
     795               &          ) / MAX( epsi10, za_s_fra(ji)  * zkappa_s(ji,0) * zg1s + (1._wp - za_s_fra(ji)) * zkappa_i(ji,0) * zg1 ) 
    811796            t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp )  ! cap t_su 
    812797         END DO 
     
    866851      !-------------------------------------------------------------------- 
    867852      ! effective conductivity and 1st layer temperature (needed by Met Office) 
     853      ! this is a conductivity at mid-layer, hence the factor 2 
    868854      DO ji = 1, npti 
    869          IF( h_i_1d(ji) > 0.1_wp ) THEN 
     855         IF( h_i_1d(ji) >= zhi_ssl ) THEN   ! zkappa_i=0 when h_i<zhi_ssl (but ztcond_i/=0) 
    870856            cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 
    871857         ELSE 
    872             cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp ! cnd_ice is capped by: cond_i/0.1m 
     858            cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / zhi_ssl ! cnd_ice is capped by: cond_i/zhi_ssl 
    873859         ENDIF 
    874          t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) 
     860         t1_ice_1d(ji) = za_s_fra(ji) * t_s_1d(ji,1) + ( 1._wp - za_s_fra(ji) ) * t_i_1d(ji,1) 
    875861      END DO 
    876862      ! 
     
    886872      DO ji = 1, npti          
    887873         !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    888          zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 
    889          IF( h_s_1d(ji) >= zhs_min ) THEN !!clem change that 
    890             t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   & 
    891                &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac ) 
     874         IF( h_s_1d(ji) >= zhs_ssl ) THEN 
     875            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) ) / & 
     876               &          ( rn_cnd_s * zh_i(ji)                + ztcond_i(ji,1) * zh_s(ji)                ) 
    892877         ELSE 
    893878            t_si_1d(ji) = t_su_1d(ji) 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icevar.F90

    r12854 r12890  
    5252   !!   ice_var_itd       : convert N-cat to M-cat 
    5353   !!   ice_var_snwfra    : fraction of ice covered by snow 
     54   !!   ice_var_snwblow   : distribute snow fall between ice and ocean 
    5455   !!---------------------------------------------------------------------- 
    5556   USE dom_oce        ! ocean space and time domain 
     
    7980   PUBLIC   ice_var_itd 
    8081   PUBLIC   ice_var_snwfra 
     82   PUBLIC   ice_var_snwblow 
    8183 
    8284   INTERFACE ice_var_itd 
     
    8688   INTERFACE ice_var_snwfra 
    8789      MODULE PROCEDURE ice_var_snwfra_1d, ice_var_snwfra_2d, ice_var_snwfra_3d 
     90   END INTERFACE 
     91 
     92   INTERFACE ice_var_snwblow 
     93      MODULE PROCEDURE ice_var_snwblow_1d, ice_var_snwblow_2d 
    8894   END INTERFACE 
    8995 
     
    10091015         pt_su(:,jl) = ptmsu(:) 
    10101016         ps_i (:,jl) = psmi (:) 
    1011          ps_i (:,jl) = psmi (:)          
    10121017      END DO 
    10131018      ! 
     
    10591064      !! 
    10601065      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1 
    1061        !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
     1066      !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
    10621067      !!               
    10631068      !!               3) Expand the filling to the empty cat between jlmin and jlmax  
     
    12931298   !!------------------------------------------------------------------- 
    12941299   SUBROUTINE ice_var_snwfra_3d( ph_s, pa_s_fra ) 
    1295       !!------------------------------------------------------------------- 
    12961300      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ph_s        ! snow thickness 
    12971301      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow 
    1298       !!------------------------------------------------------------------- 
    1299       pa_s_fra(:,:,:) = 1._wp - EXP( -0.2_wp * rhos * ph_s(:,:,:) ) 
     1302      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover 
     1303         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp 
     1304         ELSEWHERE             ; pa_s_fra = 0._wp 
     1305         END WHERE 
     1306      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style) 
     1307         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s ) 
     1308      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style) 
     1309         pa_s_fra = ph_s / ( ph_s + 0.02_wp ) 
     1310      ENDIF 
    13001311   END SUBROUTINE ice_var_snwfra_3d 
    13011312 
    13021313   SUBROUTINE ice_var_snwfra_2d( ph_s, pa_s_fra ) 
    1303       !!------------------------------------------------------------------- 
    13041314      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ph_s        ! snow thickness 
    13051315      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow 
    1306       !!------------------------------------------------------------------- 
    1307       pa_s_fra(:,:) = 1._wp - EXP( -0.2_wp * rhos * ph_s(:,:) ) 
     1316      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover 
     1317         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp 
     1318         ELSEWHERE             ; pa_s_fra = 0._wp 
     1319         END WHERE 
     1320      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style) 
     1321         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s ) 
     1322      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style) 
     1323         pa_s_fra = ph_s / ( ph_s + 0.02_wp ) 
     1324      ENDIF 
    13081325   END SUBROUTINE ice_var_snwfra_2d 
    13091326 
    13101327   SUBROUTINE ice_var_snwfra_1d( ph_s, pa_s_fra ) 
    1311       !!------------------------------------------------------------------- 
    13121328      REAL(wp), DIMENSION(:), INTENT(in   ) ::   ph_s        ! snow thickness 
    13131329      REAL(wp), DIMENSION(:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow 
    1314       !!------------------------------------------------------------------- 
    1315       pa_s_fra(:) = 1._wp - EXP( -0.2_wp * rhos * ph_s(:) ) 
     1330      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover 
     1331         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp 
     1332         ELSEWHERE             ; pa_s_fra = 0._wp 
     1333         END WHERE 
     1334      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style) 
     1335         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s ) 
     1336      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style) 
     1337         pa_s_fra = ph_s / ( ph_s + 0.02_wp ) 
     1338      ENDIF 
    13161339   END SUBROUTINE ice_var_snwfra_1d 
    13171340    
     1341   !!-------------------------------------------------------------------------- 
     1342   !! INTERFACE ice_var_snwblow 
     1343   !! 
     1344   !! ** Purpose :   Compute distribution of precip over the ice 
     1345   !! 
     1346   !!                Snow accumulation in one thermodynamic time step 
     1347   !!                snowfall is partitionned between leads and ice. 
     1348   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads 
     1349   !!                but because of the winds, more snow falls on leads than on sea ice 
     1350   !!                and a greater fraction (1-at_i)^beta of the total mass of snow  
     1351   !!                (beta < 1) falls in leads. 
     1352   !!                In reality, beta depends on wind speed,  
     1353   !!                and should decrease with increasing wind speed but here, it is  
     1354   !!                considered as a constant. an average value is 0.66 
     1355   !!-------------------------------------------------------------------------- 
     1356!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE.... 
     1357   SUBROUTINE ice_var_snwblow_2d( pin, pout ) 
     1358      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b ) 
     1359      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout 
     1360      pout = ( 1._wp - ( pin )**rn_snwblow ) 
     1361   END SUBROUTINE ice_var_snwblow_2d 
     1362 
     1363   SUBROUTINE ice_var_snwblow_1d( pin, pout ) 
     1364      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin 
     1365      REAL(wp), DIMENSION(:), INTENT(inout) :: pout 
     1366      pout = ( 1._wp - ( pin )**rn_snwblow ) 
     1367   END SUBROUTINE ice_var_snwblow_1d 
    13181368 
    13191369#else 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk.F90

    r12811 r12890  
    4747#if defined key_si3 
    4848   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif 
    49    USE icethd_dh      ! for CALL ice_thd_snwblow 
     49   USE icevar         ! for CALL ice_var_snwblow 
    5050#endif 
    5151   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
     
    890890      ! --- evaporation minus precipitation --- ! 
    891891      zsnw(:,:) = 0._wp 
    892       CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing 
     892      CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing 
    893893      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
    894894      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbccpl.F90

    r12811 r12890  
    4141#endif 
    4242#if defined key_si3 
    43    USE icethd_dh      ! for CALL ice_thd_snwblow 
     43   USE icevar         ! for CALL ice_var_snwblow 
    4444#endif 
    4545   ! 
     
    17301730 
    17311731      ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing) 
    1732       zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw ) 
     1732      zsnw(:,:) = 0._wp   ;   CALL ice_var_snwblow( ziceld, zsnw ) 
    17331733       
    17341734      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- ! 
Note: See TracChangeset for help on using the changeset viewer.