Ignore:
Timestamp:
2017-10-03T17:09:40+02:00 (3 years ago)
Author:
clem
Message:

add option Lupkes2015 for ice-atm drag

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r8422 r8585  
    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 
     
    9393   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation 
    9494   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant 
    95    REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice 
     95   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice 
    9696   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant 
    9797   ! 
     
    108108   REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements 
    109109   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements 
    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   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau) 
     114   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens) 
     115   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat) 
     116   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme) 
     117   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme) 
     118   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 
    113119 
    114120   INTEGER  ::   nblk           ! choice of the bulk algorithm 
     
    132138      !!             ***  ROUTINE sbc_blk_alloc *** 
    133139      !!------------------------------------------------------------------- 
    134       ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 
     140      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 
     141         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 
    135142      ! 
    136143      IF( lk_mpp             )   CALL mpp_sum ( sbc_blk_alloc ) 
     
    164171         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    165172         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
    166          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12 
     173         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 
    167174      !!--------------------------------------------------------------------- 
    168175      ! 
     
    255262         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac 
    256263         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))' 
     264         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12 
     265         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15 
    257266         ! 
    258267         WRITE(numout,*) 
     
    361370      REAL(wp), DIMENSION(:,:), POINTER ::   zqlw, zqsb        ! long wave and sensible heat fluxes 
    362371      REAL(wp), DIMENSION(:,:), POINTER ::   zqla, zevap       ! latent heat fluxes and evaporation 
    363       REAL(wp), DIMENSION(:,:), POINTER ::   Cd                ! transfer coefficient for momentum      (tau) 
    364       REAL(wp), DIMENSION(:,:), POINTER ::   Ch                ! transfer coefficient for sensible heat (Q_sens) 
    365       REAL(wp), DIMENSION(:,:), POINTER ::   Ce                ! tansfert coefficient for evaporation   (Q_lat) 
    366372      REAL(wp), DIMENSION(:,:), POINTER ::   zst               ! surface temperature in Kelvin 
    367       REAL(wp), DIMENSION(:,:), POINTER ::   zt_zu             ! air temperature at wind speed height 
    368       REAL(wp), DIMENSION(:,:), POINTER ::   zq_zu             ! air spec. hum.  at wind speed height 
    369373      REAL(wp), DIMENSION(:,:), POINTER ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    370374      REAL(wp), DIMENSION(:,:), POINTER ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
     
    375379      ! 
    376380      CALL wrk_alloc( jpi,jpj,   zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap ) 
    377       CALL wrk_alloc( jpi,jpj,   Cd, Ch, Ce, zst, zt_zu, zq_zu ) 
    378       CALL wrk_alloc( jpi,jpj,   zU_zu, ztpot, zrhoa ) 
     381      CALL wrk_alloc( jpi,jpj,   zst, zU_zu, ztpot, zrhoa ) 
    379382      ! 
    380383 
     
    423426      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    424427 
    425  
    426  
    427428      ! ----------------------------------------------------------------------------- ! 
    428429      !     II    Turbulent FLUXES                                                    ! 
     
    440441      ! 
    441442      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2 
    442          &                                               Cd, Ch, Ce, 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 ) 
    443444      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0 
    444          &                                               Cd, Ch, Ce, 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 ) 
    445446      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5 
    446          &                                               Cd, Ch, Ce, 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 ) 
    447448      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF 
    448          &                                               Cd, Ch, Ce, zt_zu, zq_zu, zU_zu ) 
     449         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    449450      CASE DEFAULT 
    450451         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
     
    453454      !                          ! Compute true air density : 
    454455      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 
    455          zrhoa(:,:) = rho_air( zt_zu(:,:)             , zq_zu(:,:)             , sf(jp_slp)%fnow(:,:,1) ) 
     456         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) ) 
    456457      ELSE                                      ! At zt: 
    457458         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    458459      END IF 
    459460 
    460       Cd_oce(:,:) = Cd(:,:)  ! record value of pure ocean-atm. drag (clem) 
     461!!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
     462!!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef. 
    461463 
    462464      DO jj = 1, jpj             ! tau module, i and j component 
    463465         DO ji = 1, jpi 
    464             zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd(ji,jj)   ! using bulk wind speed 
     466            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed 
    465467            taum  (ji,jj) = zztmp * wndm  (ji,jj) 
    466468            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
     
    497499      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    498500         !! q_air and t_air are given at 10m (wind reference height) 
    499          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
    500          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
     501         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
     502         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
    501503      ELSE 
    502504         !! q_air and t_air are not given at 10m (wind reference height) 
    503505         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
    504          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed 
    505          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - zt_zu(:,:) )   ! Sensible Heat ! using bulk wind speed 
     506         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 
     507         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
    506508      ENDIF 
    507509 
     
    510512 
    511513      IF(ln_ctl) THEN 
    512          CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce , clinfo2=' Ce  : ' ) 
    513          CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch , clinfo2=' Ch  : ' ) 
     514         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' ) 
     515         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' ) 
    514516         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
    515517         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' ) 
     
    563565      ! 
    564566      CALL wrk_dealloc( jpi,jpj,   zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap ) 
    565       CALL wrk_dealloc( jpi,jpj,   Cd, Ch, Ce, zst, zt_zu, zq_zu ) 
    566       CALL wrk_dealloc( jpi,jpj,   zU_zu, ztpot, zrhoa ) 
     567      CALL wrk_dealloc( jpi,jpj,   zst, zU_zu, ztpot, zrhoa ) 
    567568      ! 
    568569      IF( nn_timing == 1 )  CALL timing_stop('blk_oce') 
     
    588589      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point 
    589590      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point 
    590       REAL(wp), DIMENSION(:,:), POINTER ::   Cd               ! transfer coefficient for momentum      (tau) 
    591591      !!--------------------------------------------------------------------- 
    592592      ! 
     
    594594      ! 
    595595      CALL wrk_alloc( jpi,jpj, zrhoa ) 
    596       CALL wrk_alloc( jpi,jpj, Cd ) 
    597  
    598       Cd(:,:) = Cd_ice 
    599  
    600       ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 
    601       IF( ln_Cd_L12 ) THEN 
    602          CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 
    603       ENDIF 
    604  
    605       ! local scalars ( place there for vector optimisation purposes) 
    606       ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    607       !! 
    608       zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
    609  
    610       !!gm brutal.... 
    611       utau_ice  (:,:) = 0._wp 
    612       vtau_ice  (:,:) = 0._wp 
    613       wndm_ice  (:,:) = 0._wp 
    614       !!gm end 
    615  
    616       ! ----------------------------------------------------------------------------- ! 
    617       !    Wind components and module relative to the moving ocean ( U10m - U_ice )   ! 
    618       ! ----------------------------------------------------------------------------- ! 
     596 
     597      ! set transfer coefficients to default sea-ice values 
     598      Cd_atm(:,:) = Cd_ice 
     599      Ch_atm(:,:) = Cd_ice 
     600      Ce_atm(:,:) = Cd_ice 
     601 
     602      wndm_ice(:,:) = 0._wp      !!gm brutal.... 
     603 
     604      ! ------------------------------------------------------------ ! 
     605      !    Wind module relative to the moving ice ( U10m - U_ice )   ! 
     606      ! ------------------------------------------------------------ ! 
    619607      SELECT CASE( cp_ice_msh ) 
    620608      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
     
    622610         DO jj = 2, jpjm1 
    623611            DO ji = 2, jpim1   ! B grid : NO vector opt 
    624                ! ... scalar wind at I-point (fld being at T-point) 
    625                zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
    626                   &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
    627                zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    628                   &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    629                zwnorm_f = zrhoa(ji,jj) * Cd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    630                ! ... ice stress at I-point 
    631                utau_ice(ji,jj) = zwnorm_f * zwndi_f 
    632                vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
    633612               ! ... scalar wind at T-point (fld being at T-point) 
    634613               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
     
    639618            END DO 
    640619         END DO 
    641          CALL lbc_lnk( utau_ice, 'I', -1. ) 
    642          CALL lbc_lnk( vtau_ice, 'I', -1. ) 
    643620         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    644621         ! 
    645622      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    646          DO jj = 2, jpj 
    647             DO ji = fs_2, jpi   ! vect. opt. 
     623         DO jj = 2, jpjm1 
     624            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    648625               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
    649626               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
     
    651628            END DO 
    652629         END DO 
     630         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
     631         ! 
     632      END SELECT 
     633 
     634      ! Make ice-atm. drag dependent on ice concentration 
     635      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations 
     636         CALL Cdn10_Lupkes2012( Cd_atm ) 
     637         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical 
     638      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations 
     639         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )  
     640      ENDIF 
     641 
     642!!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef. 
     643!!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef. 
     644 
     645      ! local scalars ( place there for vector optimisation purposes) 
     646      ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
     647      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
     648 
     649      !!gm brutal.... 
     650      utau_ice  (:,:) = 0._wp 
     651      vtau_ice  (:,:) = 0._wp 
     652      !!gm end 
     653 
     654      ! ------------------------------------------------------------ ! 
     655      !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     656      ! ------------------------------------------------------------ ! 
     657      SELECT CASE( cp_ice_msh ) 
     658      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
     659         DO jj = 2, jpjm1 
     660            DO ji = 2, jpim1   ! B grid : NO vector opt 
     661               ! ... scalar wind at I-point (fld being at T-point) 
     662               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
     663                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
     664               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
     665                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
     666               ! ... ice stress at I-point 
     667               zwnorm_f = zrhoa(ji,jj) * Cd_atm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
     668               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
     669               vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
     670            END DO 
     671         END DO 
     672         CALL lbc_lnk( utau_ice, 'I', -1. ) 
     673         CALL lbc_lnk( vtau_ice, 'I', -1. ) 
     674         ! 
     675      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    653676         DO jj = 2, jpjm1 
    654677            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    655                utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     678               utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            & 
    656679                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    657                vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     680               vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            & 
    658681                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    659682            END DO 
     
    661684         CALL lbc_lnk( utau_ice, 'U', -1. ) 
    662685         CALL lbc_lnk( vtau_ice, 'V', -1. ) 
    663          CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    664686         ! 
    665687      END SELECT 
     
    700722      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3) 
    701723      REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa 
    702       REAL(wp), DIMENSION(:,:)  , POINTER ::   Cd            ! transfer coefficient for momentum      (tau) 
    703724      !!--------------------------------------------------------------------- 
    704725      ! 
     
    706727      ! 
    707728      CALL wrk_alloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    708       CALL wrk_alloc( jpi,jpj,       zrhoa) 
    709       CALL wrk_alloc( jpi,jpj, Cd ) 
    710  
    711       Cd(:,:) = Cd_ice 
    712  
    713       ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al.  2012) (clem) 
    714       IF( ln_Cd_L12 ) THEN 
    715          CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 
    716       ENDIF 
    717  
    718       ! 
    719       ! local scalars ( place there for vector optimisation purposes) 
     729      CALL wrk_alloc( jpi,jpj,       zrhoa ) 
     730      ! 
     731      ! local scalars 
    720732      zcoef_dqlw   = 4.0 * 0.95 * Stef 
    721733      zcoef_dqla   = -Ls * 11637800. * (-5897.8) 
     
    745757               ! ----------------------------! 
    746758 
    747                ! ... turbulent heat fluxes 
     759               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 
    748760               ! Sensible Heat 
    749                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) ) 
     761               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)) 
    750762               ! Latent Heat 
    751                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Cd(ji,jj) * wndm_ice(ji,jj)   & 
    752                   &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
     763               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
     764                  &                ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
    753765               ! Latent heat sensitivity for ice (Dqla/Dt) 
    754766               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
    755                   dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Cd(ji,jj) * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 
     767                  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)) 
    756768               ELSE 
    757769                  dqla_ice(ji,jj,jl) = 0._wp 
     
    759771 
    760772               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    761                z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Cd(ji,jj) * wndm_ice(ji,jj) 
     773               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) 
    762774 
    763775               ! ----------------------------! 
     
    838850      CALL wrk_dealloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    839851      CALL wrk_dealloc( jpi,jpj,       zrhoa ) 
    840       CALL wrk_dealloc( jpi,jpj, Cd ) 
    841852      ! 
    842853      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_flx') 
     
    965976 
    966977#if defined key_lim3 
     978 
    967979   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
    968980      !!---------------------------------------------------------------------- 
     
    10141026       
    10151027   END SUBROUTINE Cdn10_Lupkes2012 
     1028 
     1029 
     1030   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch ) 
     1031      !!---------------------------------------------------------------------- 
     1032      !!                      ***  ROUTINE  Cdn10_Lupkes2015  *** 
     1033      !! 
     1034      !! ** pUrpose :    1lternative turbulent transfert coefficients formulation 
     1035      !!                 between sea-ice and atmosphere with distinct momentum  
     1036      !!                 and heat coefficients depending on sea-ice concentration  
     1037      !!                 and atmospheric stability (no meltponds effect for now). 
     1038      !!                 
     1039      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015) 
     1040      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 
     1041      !!                 it considers specific skin and form drags (Andreas et al. 2010) 
     1042      !!                 to compute neutral transfert coefficients for both heat and  
     1043      !!                 momemtum fluxes. Atmospheric stability effect on transfert 
     1044      !!                 coefficient is also taken into account following Louis (1979). 
     1045      !! 
     1046      !! ** References : Lupkes et al. JGR 2015 (theory) 
     1047      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation) 
     1048      !! 
     1049      !!---------------------------------------------------------------------- 
     1050      ! 
     1051      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
     1052      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch 
     1053      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat 
     1054      ! 
     1055      ! ECHAM6 constants 
     1056      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m] 
     1057      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m] 
     1058      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m] 
     1059      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41 
     1060      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41 
     1061      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13 
     1062      REAL(wp), PARAMETER ::   zc2          = zc * zc 
     1063      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14 
     1064      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30 
     1065      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51 
     1066      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56 
     1067      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26 
     1068      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26 
     1069      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma 
     1070      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp 
     1071      ! 
     1072      INTEGER  ::   ji, jj         ! dummy loop indices 
     1073      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu 
     1074      REAL(wp) ::   zrib_o, zrib_i 
     1075      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice 
     1076      REAL(wp) ::   zChn_skin_ice, zChn_form_ice 
     1077      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw 
     1078      !!---------------------------------------------------------------------- 
     1079      ! 
     1080      ! Atmospheric and Surface Variables 
     1081      zst(:,:)     = sst_m(:,:) + rt0                                       ! convert SST from Celcius to Kelvin 
     1082      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)  , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg] 
     1083      zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] 
     1084      ! 
     1085      DO jj = 2, jpjm1 
     1086         DO ji = fs_2, fs_jpim1 
     1087            ! Virtual potential temperature [K] 
     1088            zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean 
     1089            zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice 
     1090            zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu 
     1091             
     1092            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 
     1093            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean 
     1094            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice 
     1095             
     1096            ! Momentum Neutral Transfert Coefficients (should be a constant) 
     1097            zCdn_skin_ice =         ( vkarmn                              / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2      ! Eq. 7 
     1098            zCdn_form_ice = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   &  ! Eq. 40 
     1099               &                  * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta 
     1100            zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details) 
     1101            !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32) 
     1102 
     1103            ! Heat Neutral Transfert Coefficients 
     1104            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) 
     1105            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53  
     1106                        
     1107            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 
     1108            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 
     1109            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details) 
     1110            IF( zrib_o <= 0._wp ) THEN 
     1111               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 
     1112               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * wndm(ji,jj) )   &                ! Eq. 26 
     1113                  &             )**zgamma )**z1_gamma 
     1114            ELSE 
     1115               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12 
     1116               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28 
     1117            ENDIF 
     1118             
     1119            IF( zrib_i <= 0._wp ) THEN 
     1120               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9 
     1121               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25 
     1122            ELSE 
     1123               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11 
     1124               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27 
     1125            ENDIF 
     1126             
     1127            ! Momentum Transfert Coefficients (Eq. 38) 
     1128            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  & 
     1129               &        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) ) 
     1130             
     1131            ! Heat Transfert Coefficients (Eq. 49) 
     1132            Ch(ji,jj) = zChn_skin_ice *   zfhi +  & 
     1133               &        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) ) 
     1134            ! 
     1135         END DO 
     1136      END DO 
     1137      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. ) 
     1138      ! 
     1139   END SUBROUTINE Cdn10_Lupkes2015 
     1140 
    10161141#endif 
    1017     
    10181142 
    10191143   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.