New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 8752 for branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90 – NEMO

Ignore:
Timestamp:
2017-11-20T13:54:32+01:00 (6 years ago)
Author:
dancopsey
Message:

Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) from revision 8587 to 8726.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r8738 r8752  
    1717   !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 
    1818   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine 
     19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle)  
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    4243   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 
    4344   USE icethd_dh      ! for CALL ice_thd_snwblow 
     45   USE icethd_zdf, ONLY:   rn_cnd_s  ! for blk_ice_qcn 
    4446#endif 
    4547   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
     
    6466   PUBLIC   blk_ice_tau   ! routine called in icestp module 
    6567   PUBLIC   blk_ice_flx   ! routine called in icestp module 
    66 #endif 
     68   PUBLIC   blk_ice_qcn   ! routine called in icestp module 
     69#endif  
    6770 
    6871!!Lolo: should ultimately be moved in the module with all physical constants ? 
     
    697700 
    698701 
    699    SUBROUTINE blk_ice_flx( ptsu, palb ) 
     702   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 
    700703      !!--------------------------------------------------------------------- 
    701704      !!                     ***  ROUTINE blk_ice_flx  *** 
     
    710713      !!--------------------------------------------------------------------- 
    711714      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature 
     715      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness 
     716      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness 
    712717      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies) 
    713718      !! 
     
    716721      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    717722      REAL(wp) ::   zztmp, z1_lsub           !   -      - 
     723      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice 
     724      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables 
    718725      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw         ! long wave heat flux over ice 
    719726      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb         ! sensible  heat flux over ice 
     
    722729      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3) 
    723730      REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa 
     731 
    724732      !!--------------------------------------------------------------------- 
    725733      ! 
     
    830838      CALL wrk_dealloc( jpi,jpj,   zevap, zsnw ) 
    831839 
    832       !-------------------------------------------------------------------- 
    833       ! FRACTIONs of net shortwave radiation which is not absorbed in the 
    834       ! thin surface layer and penetrates inside the ice cover 
    835       ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    836       ! 
    837       fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    838       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
    839       ! 
    840       ! 
     840      ! --- absorbed and transmitted shortwave radiation (W/m2) --- ! 
     841      ! 
     842      ! former coding was 
     843      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     844      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     845 
     846      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- ! 
     847      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value 
     848      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1 
     849 
     850      qsr_ice_tr(:,:,:) = 0._wp 
     851 
     852      DO jl = 1, jpl 
     853         DO jj = 1, jpj 
     854            DO ji = 1, jpi 
     855 
     856               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1 
     857               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 
     858 
     859               IF ( phs(ji,jj,jl) <= 0.0_wp ) THEN  
     860                  zfrqsr_tr_i  = zfr1            + zfac * zfr2 
     861               ELSE  
     862                  zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque 
     863               ENDIF 
     864 
     865               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation  
     866 
     867            END DO 
     868         END DO 
     869      END DO 
     870                                                                                       
     871 
    841872      IF(ln_ctl) THEN 
    842873         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl) 
     
    854885 
    855886   END SUBROUTINE blk_ice_flx 
     887    
     888    
     889 
     890   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi ) 
     891    
     892      !!--------------------------------------------------------------------- 
     893      !!                     ***  ROUTINE blk_ice_qcn  *** 
     894      !! 
     895      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux 
     896      !!                to force sea ice / snow thermodynamics 
     897      !!                in the case JULES coupler is emulated 
     898      !!                 
     899      !! ** Method  :   compute surface energy balance assuming neglecting heat storage 
     900      !!                following the 0-layer Semtner (1976) approach 
     901      !! 
     902      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K) 
     903      !!              - qcn_ice : surface inner conduction flux (W/m2) 
     904      !! 
     905      !!--------------------------------------------------------------------- 
     906      !! 
     907      INTEGER, INTENT(in)                        ::   k_monocat                 ! single-category option 
     908 
     909      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   ptsu                      ! sea ice / snow surface temperature 
     910 
     911      REAL(wp), DIMENSION(:,:)  , INTENT(in)     ::   ptb                       ! sea ice base temperature 
     912      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phs                       ! snow thickness 
     913      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phi                       ! sea ice thickness 
     914       
     915      INTEGER                                    ::   ji, jj, jl                ! dummy loop indices 
     916      INTEGER                                    ::   iter                      ! 
     917      REAL(wp)                                   ::   zfac, zfac2, zfac3        ! dummy factors 
     918      REAL(wp)                                   ::   zkeff_h, ztsu             ! 
     919      REAL(wp)                                   ::   zqc, zqnet                ! 
     920      REAL(wp)                                   ::   zhe, zqa0                 ! 
     921            
     922      INTEGER , PARAMETER                        ::   nit = 10                  ! number of iterations 
     923      REAL(wp), PARAMETER                        ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction 
     924 
     925      REAL(wp), DIMENSION(:,:,:), POINTER        ::   zgfac                    ! enhanced conduction factor 
     926       
     927      !!--------------------------------------------------------------------- 
     928 
     929      IF( nn_timing == 1 )  CALL timing_start('blk_ice_qcn') 
     930      ! 
     931      CALL wrk_alloc( jpi,jpj,jpl,   zgfac ) 
     932 
     933      ! -------------------------------------! 
     934      !      I   Enhanced conduction factor  ! 
     935      ! -------------------------------------! 
     936      ! 
     937      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3) 
     938      ! Fichefet and Morales Maqueda, JGR 1997 
     939      ! 
     940      zgfac(:,:,:) = 1._wp 
     941       
     942      SELECT CASE ( k_monocat )  
     943       
     944         CASE ( 1 , 3 ) 
     945       
     946            zfac    =   1._wp /  ( rn_cnd_s + rcdic ) 
     947            zfac2   =   EXP(1._wp) * 0.5_wp * zepsilon 
     948            zfac3   =   2._wp / zepsilon 
     949       
     950            DO jl = 1, jpl                 
     951               DO jj = 1 , jpj 
     952                  DO ji = 1, jpi 
     953                     !                                ! Effective thickness 
     954                     zhe               =   ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac 
     955                     !                                ! Enhanced conduction factor 
     956                     IF( zhe >=  zfac2 )  & 
     957                        zgfac(ji,jj,jl)   =   MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) ) 
     958                  END DO 
     959               END DO 
     960            END DO 
     961       
     962      END SELECT 
     963       
     964      ! -------------------------------------------------------------! 
     965      !      II   Surface temperature and conduction flux            ! 
     966      ! -------------------------------------------------------------! 
     967 
     968      zfac   =   rcdic * rn_cnd_s  
     969                                      ! ========================== ! 
     970      DO jl = 1, jpl                  !  Loop over ice categories  ! 
     971         !                            ! ========================== ! 
     972         DO jj = 1 , jpj 
     973            DO ji = 1, jpi 
     974               !                      ! Effective conductivity of the snow-ice system divided by thickness 
     975               zkeff_h   =   zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) ) 
     976                                      ! Store initial surface temperature 
     977               ztsu      =   ptsu(ji,jj,jl) 
     978                                      ! Net initial atmospheric heat flux 
     979               zqa0      =   qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl) 
     980                
     981               DO iter = 1, nit       ! --- Iteration loop 
     982                
     983                   !                  ! Conduction heat flux through snow-ice system (>0 downwards) 
     984                   zqc   =   zkeff_h * ( ztsu - ptb(ji,jj) ) 
     985                   !                  ! Surface energy budget 
     986                   zqnet =   zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc 
     987                   !                  ! Temperature update 
     988                   ztsu  =   ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) 
     989                
     990               END DO 
     991                
     992               ptsu(ji,jj,jl)    = MIN( rt0, ztsu ) 
     993                
     994               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 
     995                
     996            END DO 
     997         END DO 
     998          
     999      END DO  
     1000       
     1001      CALL wrk_dealloc( jpi,jpj,jpl,   zgfac ) 
     1002       
     1003      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_qcn') 
     1004 
     1005   END SUBROUTINE blk_ice_qcn 
    8561006    
    8571007#endif 
     
    10931243      zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] 
    10941244      ! 
    1095 !!      DO jj = 2, jpjm1 
    1096 !!         DO ji = fs_2, fs_jpim1 
    1097       DO jj = 1, jpj 
    1098          DO ji = 1, jpi 
     1245      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility 
     1246         DO ji = fs_2, fs_jpim1 
    10991247            ! Virtual potential temperature [K] 
    11001248            zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean 
     
    11401288         END DO 
    11411289      END DO 
    1142 !!      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. ) 
     1290      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. ) 
    11431291      ! 
    11441292   END SUBROUTINE Cdn10_Lupkes2015 
Note: See TracChangeset for help on using the changeset viewer.