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 8637 for branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/SBC/sbcblk.F90 – NEMO

Ignore:
Timestamp:
2017-10-18T19:14:32+02:00 (6 years ago)
Author:
gm
Message:

#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8626

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/SBC/sbcblk.F90

    r8586 r8637  
    4040   USE lib_fortran    ! to use key_nosignedzero 
    4141#if defined key_lim3 
    42    USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b 
     42   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 
    4343   USE icethd_dh      ! for CALL ice_thd_snwblow 
    4444#endif 
     
    9292   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation 
    9393   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant 
    94    REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice 
     94   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice 
    9595   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant 
    9696   ! 
     
    108108   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements 
    109109!!gm ref namelist initialize it so remove the setting to false below 
    110    LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012) 
     110   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012) 
     111   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 
    111112   ! 
    112    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
     113!!cr   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
     114   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau) 
     115   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens) 
     116   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat) 
     117   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme) 
     118   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme) 
     119   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 
    113120 
    114121   INTEGER  ::   nblk           ! choice of the bulk algorithm 
     
    132139      !!             ***  ROUTINE sbc_blk_alloc *** 
    133140      !!------------------------------------------------------------------- 
    134       ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 
     141!!cr      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 
     142      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 
     143         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 
    135144      ! 
    136145      IF( lk_mpp             )   CALL mpp_sum ( sbc_blk_alloc ) 
     
    164173         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    165174         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
    166          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12 
     175         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 
    167176      !!--------------------------------------------------------------------- 
    168177      ! 
     
    255264         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac 
    256265         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))' 
     266         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12 
     267         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15 
    257268         ! 
    258269         WRITE(numout,*) 
     
    361372      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes 
    362373      REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation 
    363       REAL(wp), DIMENSION(jpi,jpj) ::   zCd               ! transfer coefficient for momentum      (tau) 
    364       REAL(wp), DIMENSION(jpi,jpj) ::   zCh               ! transfer coefficient for sensible heat (Q_sens) 
    365       REAL(wp), DIMENSION(jpi,jpj) ::   zCe               ! tansfert coefficient for evaporation   (Q_lat) 
    366374      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
    367       REAL(wp), DIMENSION(jpi,jpj) ::   zt_zu             ! air temperature at wind speed height 
    368       REAL(wp), DIMENSION(jpi,jpj) ::   zq_zu             ! air spec. hum.  at wind speed height 
    369375      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    370376      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
     
    418424      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    419425 
    420  
    421  
    422426      ! ----------------------------------------------------------------------------- ! 
    423427      !     II    Turbulent FLUXES                                                    ! 
     
    435439      ! 
    436440      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2 
    437          &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 
     441         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    438442      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0 
    439          &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 
     443         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    440444      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5 
    441          &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 
     445         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    442446      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF 
    443          &                                              zCd, zCh, zCe, zt_zu, zq_zu, zU_zu ) 
     447         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    444448      CASE DEFAULT 
    445449         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
     
    448452      !                          ! Compute true air density : 
    449453      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 
    450          zrhoa(:,:) = rho_air( zt_zu(:,:)             , zq_zu(:,:)             , sf(jp_slp)%fnow(:,:,1) ) 
     454         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) ) 
    451455      ELSE                                      ! At zt: 
    452456         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    453457      END IF 
    454458 
    455       Cd_oce(:,:) = zCd(:,:)     ! record value of pure ocean-atm. drag (clem) 
     459!!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
     460!!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef. 
    456461 
    457462      DO jj = 1, jpj             ! tau module, i and j component 
    458463         DO ji = 1, jpi 
    459             zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * zCd(ji,jj)   ! using bulk wind speed 
     464            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed 
    460465            taum  (ji,jj) = zztmp * wndm  (ji,jj) 
    461466            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
     
    492497      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    493498         !! q_air and t_air are given at 10m (wind reference height) 
    494          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*zCe(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) )  ! Evaporation, using bulk wind speed 
    495          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*zCh(:,:)*(zst(:,:) - ztpot(:,:)               )   ! Sensible Heat, using bulk wind speed 
     499         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
     500         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
    496501      ELSE 
    497502         !! q_air and t_air are not given at 10m (wind reference height) 
    498503         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
    499          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*zCe(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed 
    500          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*zCh(:,:)*(zst(:,:) - zt_zu(:,:) )   ! Sensible Heat ! using bulk wind speed 
     504         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 
     505         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
    501506      ENDIF 
    502507 
     
    505510 
    506511      IF(ln_ctl) THEN 
    507          CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=zCe , clinfo2=' Ce  : ' ) 
    508          CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zCh , clinfo2=' Ch  : ' ) 
     512         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' ) 
     513         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' ) 
    509514         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
    510515         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' ) 
     
    574579      !!--------------------------------------------------------------------- 
    575580      INTEGER  ::   ji, jj    ! dummy loop indices 
    576       REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f      ! relative wind module and components at F-point 
    577       REAL(wp) ::   zwndi_t , zwndj_t                ! relative wind components at T-point 
    578       REAL(wp), DIMENSION(jpi,jpj) ::   zCd, zrhoa   ! transfer coefficient for momentum      (tau) 
     581      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point 
     582      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point 
     583      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau) 
    579584      !!--------------------------------------------------------------------- 
    580585      ! 
    581586      IF( ln_timing )   CALL timing_start('blk_ice_tau') 
    582587      ! 
    583       IF( ln_Cd_L12 ) THEN 
    584          CALL Cdn10_Lupkes2012( zCd )  ! air-ice drag = F(ice concentration) (see Lupkes et al., 2012) 
    585       ELSE 
    586          zCd(:,:) = Cd_ice             ! constant air-ice drag 
    587       ENDIF 
    588  
    589       ! local scalars ( place there for vector optimisation purposes) 
    590       ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    591       !! 
    592       zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
    593  
    594       !!gm brutal.... 
    595       utau_ice  (:,:) = 0._wp 
    596       vtau_ice  (:,:) = 0._wp 
    597       wndm_ice  (:,:) = 0._wp 
    598       !!gm end 
    599  
    600       ! ----------------------------------------------------------------------------- ! 
    601       !    Wind components and module relative to the moving ocean ( U10m - U_ice )   ! 
    602       ! ----------------------------------------------------------------------------- ! 
     588      ! set transfer coefficients to default sea-ice values 
     589      Cd_atm(:,:) = Cd_ice 
     590      Ch_atm(:,:) = Cd_ice 
     591      Ce_atm(:,:) = Cd_ice 
     592 
     593      wndm_ice(:,:) = 0._wp      !!gm brutal.... 
     594 
     595      ! ------------------------------------------------------------ ! 
     596      !    Wind module relative to the moving ice ( U10m - U_ice )   ! 
     597      ! ------------------------------------------------------------ ! 
    603598      SELECT CASE( cp_ice_msh ) 
    604599      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
     
    606601         DO jj = 2, jpjm1 
    607602            DO ji = 2, jpim1   ! B grid : NO vector opt 
    608                ! ... scalar wind at I-point (fld being at T-point) 
    609                zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
    610                   &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
    611                zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    612                   &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    613                zwnorm_f = zrhoa(ji,jj) * zCd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    614                ! ... ice stress at I-point 
    615                utau_ice(ji,jj) = zwnorm_f * zwndi_f 
    616                vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
    617603               ! ... scalar wind at T-point (fld being at T-point) 
    618604               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
     
    623609            END DO 
    624610         END DO 
    625          CALL lbc_lnk( utau_ice, 'I', -1. ) 
    626          CALL lbc_lnk( vtau_ice, 'I', -1. ) 
    627611         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    628612         ! 
    629613      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    630          DO jj = 2, jpj 
    631             DO ji = fs_2, jpi   ! vect. opt. 
     614         DO jj = 2, jpjm1 
     615            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    632616               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
    633617               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
     
    635619            END DO 
    636620         END DO 
     621         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
     622         ! 
     623      END SELECT 
     624 
     625      ! Make ice-atm. drag dependent on ice concentration 
     626      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations 
     627         CALL Cdn10_Lupkes2012( Cd_atm ) 
     628         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical 
     629      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations 
     630         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )  
     631      ENDIF 
     632 
     633!!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef. 
     634!!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef. 
     635 
     636      ! local scalars ( place there for vector optimisation purposes) 
     637      ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
     638      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
     639 
     640      !!gm brutal.... 
     641      utau_ice  (:,:) = 0._wp 
     642      vtau_ice  (:,:) = 0._wp 
     643      !!gm end 
     644 
     645      ! ------------------------------------------------------------ ! 
     646      !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     647      ! ------------------------------------------------------------ ! 
     648      SELECT CASE( cp_ice_msh ) 
     649      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
     650         DO jj = 2, jpjm1 
     651            DO ji = 2, jpim1   ! B grid : NO vector opt 
     652               ! ... scalar wind at I-point (fld being at T-point) 
     653               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
     654                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
     655               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
     656                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
     657               ! ... ice stress at I-point 
     658               zwnorm_f = zrhoa(ji,jj) * Cd_atm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
     659               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
     660               vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
     661            END DO 
     662         END DO 
     663         CALL lbc_lnk( utau_ice, 'I', -1. ) 
     664         CALL lbc_lnk( vtau_ice, 'I', -1. ) 
     665         ! 
     666      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    637667         DO jj = 2, jpjm1 
    638668            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    639                utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * zCd(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     669               utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            & 
    640670                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    641                vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * zCd(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     671               vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            & 
    642672                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    643673            END DO 
     
    645675         CALL lbc_lnk( utau_ice, 'U', -1. ) 
    646676         CALL lbc_lnk( vtau_ice, 'V', -1. ) 
    647          CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    648677         ! 
    649678      END SELECT 
     
    684713      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3) 
    685714      REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa 
    686       REAL(wp), DIMENSION(jpi,jpj)     ::   zCd            ! transfer coefficient for momentum      (tau) 
    687715      !!--------------------------------------------------------------------- 
    688716      ! 
    689717      IF( ln_timing )   CALL timing_start('blk_ice_flx') 
    690718      ! 
    691       IF( ln_Cd_L12 ) THEN 
    692          CALL Cdn10_Lupkes2012( zCd )  ! air-ice drag = F(ice concentration) (see Lupkes et al., 2012) 
    693       ELSE 
    694          zCd(:,:) = Cd_ice             ! constant air-ice drag 
    695       ENDIF 
    696  
    697       ! 
    698       ! local scalars ( place there for vector optimisation purposes) 
     719      ! 
     720      ! local scalars 
    699721      zcoef_dqlw   = 4.0 * 0.95 * Stef 
    700722      zcoef_dqla   = -Ls * 11637800. * (-5897.8) 
     
    724746               ! ----------------------------! 
    725747 
    726                ! ... turbulent heat fluxes 
     748               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 
    727749               ! Sensible Heat 
    728                z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * zCd(ji,jj) * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
     750               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1)) 
    729751               ! Latent Heat 
    730                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * zCd(ji,jj) * wndm_ice(ji,jj)   & 
    731                   &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
     752               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
     753                  &                ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
    732754               ! Latent heat sensitivity for ice (Dqla/Dt) 
    733755               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
    734                   dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * zCd(ji,jj) * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 
     756                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) / zst2 * EXP(-5897.8 / ptsu(ji,jj,jl)) 
    735757               ELSE 
    736758                  dqla_ice(ji,jj,jl) = 0._wp 
     
    738760 
    739761               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    740                z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * zCd(ji,jj) * wndm_ice(ji,jj) 
     762               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) 
    741763 
    742764               ! ----------------------------! 
     
    935957#if defined key_lim3 
    936958 
    937    SUBROUTINE Cdn10_Lupkes2012( pCd ) 
     959   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
    938960      !!---------------------------------------------------------------------- 
    939961      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
     
    965987      !! 
    966988      !!---------------------------------------------------------------------- 
    967       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pCd   ! air-ice drag coefficient 
     989      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
    968990      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp 
    969991      REAL(wp), PARAMETER ::   znu   = 1._wp 
     
    973995      !!---------------------------------------------------------------------- 
    974996      zcoef = znu + 1._wp / ( 10._wp * zbeta ) 
    975       ! 
     997 
    976998      ! generic drag over a cell partly covered by ice 
    977       !!pCd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag 
    978       !!   &       Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag 
    979       !!   &       zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology 
     999      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag 
     1000      !!   &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag 
     1001      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology 
    9801002 
    9811003      ! ice-atm drag 
    982       pCd(:,:) = Cd_ice +  &                                                         ! pure ice drag 
    983          &       zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
    984       ! 
     1004      Cd(:,:) = Cd_ice +  &                                                         ! pure ice drag 
     1005         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology 
     1006       
    9851007   END SUBROUTINE Cdn10_Lupkes2012 
     1008 
     1009 
     1010   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch ) 
     1011      !!---------------------------------------------------------------------- 
     1012      !!                      ***  ROUTINE  Cdn10_Lupkes2015  *** 
     1013      !! 
     1014      !! ** pUrpose :    1lternative turbulent transfert coefficients formulation 
     1015      !!                 between sea-ice and atmosphere with distinct momentum  
     1016      !!                 and heat coefficients depending on sea-ice concentration  
     1017      !!                 and atmospheric stability (no meltponds effect for now). 
     1018      !!                 
     1019      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015) 
     1020      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 
     1021      !!                 it considers specific skin and form drags (Andreas et al. 2010) 
     1022      !!                 to compute neutral transfert coefficients for both heat and  
     1023      !!                 momemtum fluxes. Atmospheric stability effect on transfert 
     1024      !!                 coefficient is also taken into account following Louis (1979). 
     1025      !! 
     1026      !! ** References : Lupkes et al. JGR 2015 (theory) 
     1027      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation) 
     1028      !! 
     1029      !!---------------------------------------------------------------------- 
     1030      ! 
     1031      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
     1032      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch 
     1033      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat 
     1034      ! 
     1035      ! ECHAM6 constants 
     1036      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m] 
     1037      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m] 
     1038      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m] 
     1039      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41 
     1040      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41 
     1041      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13 
     1042      REAL(wp), PARAMETER ::   zc2          = zc * zc 
     1043      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14 
     1044      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30 
     1045      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51 
     1046      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56 
     1047      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26 
     1048      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26 
     1049      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma 
     1050      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp 
     1051      ! 
     1052      INTEGER  ::   ji, jj         ! dummy loop indices 
     1053      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu 
     1054      REAL(wp) ::   zrib_o, zrib_i 
     1055      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice 
     1056      REAL(wp) ::   zChn_skin_ice, zChn_form_ice 
     1057      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw 
     1058      REAL(wp) ::   zCdn_form_tmp 
     1059      !!---------------------------------------------------------------------- 
     1060 
     1061      ! Momentum Neutral Transfert Coefficients (should be a constant) 
     1062      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40 
     1063      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7 
     1064      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details) 
     1065      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32) 
     1066 
     1067      ! Heat Neutral Transfert Coefficients 
     1068      zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) )   ! Eq. 50 + Eq. 52 (cf Lupkes email for details) 
     1069      
     1070      ! Atmospheric and Surface Variables 
     1071      zst(:,:)     = sst_m(:,:) + rt0                                       ! convert SST from Celcius to Kelvin 
     1072      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)  , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg] 
     1073      zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] 
     1074      ! 
     1075      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility 
     1076         DO ji = fs_2, fs_jpim1 
     1077            ! Virtual potential temperature [K] 
     1078            zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean 
     1079            zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice 
     1080            zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu 
     1081             
     1082            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 
     1083            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean 
     1084            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice 
     1085             
     1086            ! Momentum and Heat Neutral Transfert Coefficients 
     1087            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40 
     1088            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53  
     1089                        
     1090            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 
     1091            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 
     1092            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details) 
     1093            IF( zrib_o <= 0._wp ) THEN 
     1094               zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) )  ! Eq. 10 
     1095               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26 
     1096                  &             )**zgamma )**z1_gamma 
     1097            ELSE 
     1098               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12 
     1099               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28 
     1100            ENDIF 
     1101             
     1102            IF( zrib_i <= 0._wp ) THEN 
     1103               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9 
     1104               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25 
     1105            ELSE 
     1106               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11 
     1107               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27 
     1108            ENDIF 
     1109             
     1110            ! Momentum Transfert Coefficients (Eq. 38) 
     1111            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  & 
     1112               &        zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 
     1113             
     1114            ! Heat Transfert Coefficients (Eq. 49) 
     1115            Ch(ji,jj) = zChn_skin_ice *   zfhi +  & 
     1116               &        zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 
     1117            ! 
     1118         END DO 
     1119      END DO 
     1120      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. ) 
     1121      ! 
     1122   END SUBROUTINE Cdn10_Lupkes2015 
     1123 
    9861124#endif 
    987     
    9881125 
    9891126   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.