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 7355 for branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90 – NEMO

Ignore:
Timestamp:
2016-11-28T18:21:42+01:00 (7 years ago)
Author:
flavoni
Message:

merge branch dev_CNRS_2016 & dev_CNRS_AGRIF_LIM3

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r7282 r7355  
    4040   USE lib_fortran    ! to use key_nosignedzero 
    4141#if defined key_lim3 
    42    USE ice     , ONLY :   u_ice, v_ice, jpl, pfrld, a_i_b 
     42   USE ice     , ONLY :   u_ice, v_ice, jpl, pfrld, a_i_b, at_i_b 
    4343   USE limthd_dh      ! for CALL lim_thd_snwblow 
    4444#elif defined key_lim2 
     
    9393 
    9494   !                                             !!! Bulk parameters 
    95    REAL(wp), PARAMETER ::   cpa  = 1000.5         ! specific heat of air (only used for ice fluxes now...) 
    96    REAL(wp), PARAMETER ::   Ls   =    2.839e6     ! latent heat of sublimation 
    97    REAL(wp), PARAMETER ::   Stef =    5.67e-8     ! Stefan Boltzmann constant 
    98    REAL(wp), PARAMETER ::   Cice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice 
    99    REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be constant 
     95   REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air (only used for ice fluxes now...) 
     96   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation 
     97   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant 
     98   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice 
     99   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant 
    100100   ! 
    101101   !                           !!* Namelist namsbc_blk : bulk parameters 
     
    111111   REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements 
    112112   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements 
     113   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012) 
     114   ! 
     115   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
    113116 
    114117   INTEGER  ::   nblk           ! choice of the bulk algorithm 
     
    128131CONTAINS 
    129132 
     133   INTEGER FUNCTION sbc_blk_alloc() 
     134      !!------------------------------------------------------------------- 
     135      !!             ***  ROUTINE sbc_blk_alloc *** 
     136      !!------------------------------------------------------------------- 
     137      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 
     138      ! 
     139      IF( lk_mpp             )   CALL mpp_sum ( sbc_blk_alloc ) 
     140      IF( sbc_blk_alloc /= 0 )   CALL ctl_warn('sbc_blk_alloc: failed to allocate arrays') 
     141   END FUNCTION sbc_blk_alloc 
     142 
    130143   SUBROUTINE sbc_blk_init 
    131144      !!--------------------------------------------------------------------- 
     
    153166         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                & 
    154167         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    155          &                 cn_dir , ln_taudif, rn_zqt, rn_zu, rn_pfac, rn_efac, rn_vfac 
    156       !!--------------------------------------------------------------------- 
     168         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
     169         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12 
     170      !!--------------------------------------------------------------------- 
     171      ! 
     172      !                                      ! allocate sbc_blk_core array 
     173      IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 
    157174      ! 
    158175      !                             !** read bulk namelist   
     
    425442      END IF 
    426443 
     444      Cd_oce(:,:) = Cd(:,:)  ! record value of pure ocean-atm. drag (clem) 
     445 
    427446      DO jj = 1, jpj             ! tau module, i and j component 
    428447         DO ji = 1, jpi 
     
    549568      INTEGER  ::   ji, jj    ! dummy loop indices 
    550569      ! 
    551       REAL(wp), DIMENSION(:,:)  , POINTER :: zcoef_wnorm 
     570      REAL(wp), DIMENSION(:,:)  , POINTER :: zrhoa 
    552571      ! 
    553572      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point 
    554573      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point 
     574      REAL(wp), DIMENSION(:,:), POINTER ::   Cd               ! transfer coefficient for momentum      (tau) 
    555575      !!--------------------------------------------------------------------- 
    556576      ! 
    557577      IF( nn_timing == 1 )  CALL timing_start('blk_ice_tau') 
    558578      ! 
    559       CALL wrk_alloc( jpi,jpj,   zcoef_wnorm ) 
     579      CALL wrk_alloc( jpi,jpj, zrhoa ) 
     580      CALL wrk_alloc( jpi,jpj, Cd ) 
     581 
     582      Cd(:,:) = Cd_ice 
     583 
     584      ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 
     585#if defined key_lim3 
     586      IF( ln_Cd_L12 ) THEN 
     587         CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 
     588      ENDIF 
     589#endif 
    560590 
    561591      ! local scalars ( place there for vector optimisation purposes) 
    562592      ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    563593      !! 
    564       zcoef_wnorm (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
    565       zcoef_wnorm (:,:) =  Cice * zcoef_wnorm (:,:) 
     594      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
    566595 
    567596      !!gm brutal.... 
     
    584613               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    585614                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    586                zwnorm_f = zcoef_wnorm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
     615               zwnorm_f = zrhoa(ji,jj) * Cd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    587616               ! ... ice stress at I-point 
    588617               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
     
    610639         DO jj = 2, jpjm1 
    611640            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    612                utau_ice(ji,jj) = 0.5*zcoef_wnorm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     641               utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
    613642                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    614                vtau_ice(ji,jj) = 0.5*zcoef_wnorm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     643               vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
    615644                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    616645            END DO 
     
    626655         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ') 
    627656      ENDIF 
    628  
    629       CALL wrk_dealloc( jpi,jpj,   zcoef_wnorm ) 
    630657 
    631658      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_tau') 
     
    659686      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3) 
    660687      REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa 
     688      REAL(wp), DIMENSION(:,:)  , POINTER ::   Cd            ! transfer coefficient for momentum      (tau) 
    661689      !!--------------------------------------------------------------------- 
    662690      ! 
     
    665693      CALL wrk_alloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    666694      CALL wrk_alloc( jpi,jpj,       zrhoa) 
     695      CALL wrk_alloc( jpi,jpj, Cd ) 
     696 
     697      Cd(:,:) = Cd_ice 
     698 
     699      ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al.  2012) (clem) 
     700#if defined key_lim3 
     701      IF( ln_Cd_L12 ) THEN 
     702         CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 
     703      ENDIF 
     704#endif 
     705 
    667706      ! 
    668707      ! local scalars ( place there for vector optimisation purposes) 
    669708      zcoef_dqlw   = 4.0 * 0.95 * Stef 
    670       zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
     709      zcoef_dqla   = -Ls * 11637800. * (-5897.8) 
    671710      ! 
    672711      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     
    696735               ! ... turbulent heat fluxes 
    697736               ! Sensible Heat 
    698                z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
     737               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Cd(ji,jj) * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    699738               ! Latent Heat 
    700                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Cice * wndm_ice(ji,jj)   & 
     739               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Cd(ji,jj) * wndm_ice(ji,jj)   & 
    701740                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    702741               ! Latent heat sensitivity for ice (Dqla/Dt) 
     
    708747 
    709748               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    710                z_dqsb(ji,jj,jl) = zrhoa(ji,jj)*cpa*Cice * wndm_ice(ji,jj) 
     749               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Cd(ji,jj) * wndm_ice(ji,jj) 
    711750 
    712751               ! ----------------------------! 
     
    782821      CALL wrk_dealloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    783822      CALL wrk_dealloc( jpi,jpj,       zrhoa ) 
     823      CALL wrk_dealloc( jpi,jpj, Cd ) 
    784824      ! 
    785825      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_flx') 
     
    906946   END FUNCTION L_vap 
    907947 
     948 
     949#if defined key_lim3 
     950   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
     951      !!---------------------------------------------------------------------- 
     952      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
     953      !! 
     954      !! ** Purpose :    Recompute the ice-atm drag at 10m height to make 
     955      !!                 it dependent on edges at leads, melt ponds and flows. 
     956      !!                 After some approximations, this can be resumed to a dependency 
     957      !!                 on ice concentration. 
     958      !!                 
     959      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50) 
     960      !!                 with the highest level of approximation: level4, eq.(59) 
     961      !!                 The generic drag over a cell partly covered by ice can be re-written as follows: 
     962      !! 
     963      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu 
     964      !! 
     965      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59) 
     966      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59) 
     967      !!                    A is the concentration of ice minus melt ponds (if any) 
     968      !! 
     969      !!                 This new drag has a parabolic shape (as a function of A) starting at 
     970      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5  
     971      !!                 and going down to Cdi(say 1.4e-3) for A=1 
     972      !! 
     973      !!                 It is theoretically applicable to all ice conditions !(not only MIZ) 
     974      !!                 => see Lupkes et al (2013) 
     975      !! 
     976      !! ** References : Lupkes et al. JGR 2012 (theory) 
     977      !!                 Lupkes et al. GRL 2013 (application to GCM) 
     978      !! 
     979      !!---------------------------------------------------------------------- 
     980      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
     981      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp 
     982      REAL(wp), PARAMETER ::   znu   = 1._wp 
     983      REAL(wp), PARAMETER ::   zmu   = 1._wp 
     984      REAL(wp), PARAMETER ::   zbeta = 1._wp 
     985      REAL(wp)            ::   zcoef 
     986      !!---------------------------------------------------------------------- 
     987      zcoef = znu + 1._wp / ( 10._wp * zbeta ) 
     988 
     989      ! generic drag over a cell partly covered by ice 
     990      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  & !! pure ocean drag 
     991      !!   &      Cd_ice      *           at_i_b(:,:)   +  & !! pure ice drag 
     992      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * !at_i_b(:,:)**zmu   ! change due to sea-ice morphology 
     993 
     994      ! ice-atm drag 
     995      Cd(:,:) = Cd_ice +  & ! pure ice drag 
     996         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)   ! change due to sea-ice morphology 
     997       
     998   END SUBROUTINE Cdn10_Lupkes2012 
     999#endif 
     1000    
     1001 
    9081002   !!====================================================================== 
    9091003END MODULE sbcblk 
Note: See TracChangeset for help on using the changeset viewer.