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 8585 – NEMO

Changeset 8585


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

add option Lupkes2015 for ice-atm drag

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r8565 r8585  
    130130&namthd_zdf     !   Ice heat diffusion 
    131131!------------------------------------------------------------------------------ 
    132    ln_zdf_Beer      = .true.          !  Heat diffusion follows a Beer law 
     132   ln_zdf_BL99      = .true.          !  Heat diffusion follows Bitz and Lipscomb 1999 
    133133   ln_cndi_U64      = .false.         !  sea ice thermal conductivity: k = k0 + beta.S/T            (Untersteiner, 1964) 
    134134   ln_cndi_P07      = .true.          !  sea ice thermal conductivity: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007) 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ref

    r8563 r8585  
    258258   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity 
    259259                           !  in the calculation of the wind stress (0.=absolute winds or 1.=relative winds) 
    260    ln_Cd_L12   = .false.   !  Modify the drag ice-atm and oce-atm depending on ice concentration 
    261                            !  This parameterization is from Lupkes et al. (JGR 2012) 
     260   ln_Cd_L12   = .false.   !  Modify the drag ice-atm depending on ice concentration with Lupkes 2012 
     261   ln_Cd_L15   = .false.   !  Modify the drag ice-atm depending on ice concentration with Lupkes 2015 
    262262/ 
    263263!----------------------------------------------------------------------- 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90

    r8565 r8585  
    3131 
    3232   !!** namelist (namthd_zdf) ** 
    33    LOGICAL  ::   ln_zdf_Beer      ! Heat diffusion follows a Beer Law 
     33   LOGICAL  ::   ln_zdf_BL99      ! Heat diffusion follows Bitz and Lipscomb (1999) 
    3434   LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964) 
    3535   LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007) 
    … …  
    677677      INTEGER  ::   ios   ! Local integer output status for namelist read 
    678678      !! 
    679       NAMELIST/namthd_zdf/ ln_zdf_Beer, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i, ln_dqns_i 
     679      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i, ln_dqns_i 
    680680      !!------------------------------------------------------------------- 
    681681      ! 
    … …  
    694694         WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
    695695         WRITE(numout,*) '   Namelist namthd_zdf:' 
    696          WRITE(numout,*) '      Diffusion follows a Beer Law                            ln_zdf_Beer  = ', ln_zdf_Beer 
     696         WRITE(numout,*) '      Diffusion follows a Bitz and Lipscomb (1999)            ln_zdf_BL99  = ', ln_zdf_BL99 
    697697         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64 
    698698         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07 

  • 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   !!====================================================================== 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90

    r7646 r8585  
    5252   !! COARE own values for given constants: 
    5353   REAL(wp), PARAMETER :: & 
    54       &   zi0     = 600.,  &   !: scale height of the atmospheric boundary layer...1 
     54      &   zi0     = 600., &   !: scale height of the atmospheric boundary layer...1 
    5555      &  Beta0    = 1.25, &   !: gustiness parameter 
    5656      &  rctv0    = 0.608     !: constant to obtain virtual temperature... 
    … …  
    6060 
    6161   SUBROUTINE turb_coare( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 
    62       &                   Cd, Ch, Ce, t_zu, q_zu, U_blk ) 
     62      &                   Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
     63      &                   Cdn, Chn, Cen                       ) 
     64 
    6365      !!---------------------------------------------------------------------- 
    6466      !!                      ***  ROUTINE  turb_coare  *** 
    … …  
    106108      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    107109      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
     110      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    108111      ! 
    109112      INTEGER :: j_itt 
    … …  
    246249      Ce   = ztmp0*q_star/dq_zu 
    247250      ! 
     251      ztmp1 = zu + z0 
     252      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 
     253      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 
     254      Cen = Chn 
     255      ! 
    248256      CALL wrk_dealloc( jpi,jpj, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu ) 
    249257      CALL wrk_dealloc( jpi,jpj, znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 ) 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare3p5.F90

    r7646 r8585  
    5858CONTAINS 
    5959 
    60    SUBROUTINE turb_coare3p5( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 
    61       &                   Cd, Ch, Ce, t_zu, q_zu, U_blk ) 
     60   SUBROUTINE turb_coare3p5( zt, zu, sst, t_zt, ssq, q_zt, U_zu,  & 
     61      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,       & 
     62      &                      Cdn, Chn, Cen                        ) 
     63 
    6264      !!---------------------------------------------------------------------------------- 
    6365      !!                      ***  ROUTINE  turb_coare3p5  *** 
    … …  
    105107      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu             [kg/kg] 
    106108      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
     109      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    107110      ! 
    108111      INTEGER :: j_itt 
    … …  
    252255      Ce   = ztmp0*q_star/dq_zu 
    253256      ! 
     257      ztmp1 = zu + z0 
     258      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 
     259      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 
     260      Cen = Chn 
     261      ! 
    254262      CALL wrk_dealloc( jpi,jpj, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu ) 
    255263      CALL wrk_dealloc( jpi,jpj, znu_a, z0, z0t, ztmp0, ztmp1, ztmp2 ) 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ecmwf.F90

    r7646 r8585  
    6464 
    6565   SUBROUTINE TURB_ECMWF( zt, zu, sst, t_zt, ssq , q_zt , U_zu,   & 
    66       &                   Cd, Ch, Ce , t_zu, q_zu, U_blk ) 
     66      &                   Cd, Ch, Ce , t_zu, q_zu, U_blk,         & 
     67      &                   Cdn, Chn, Cen                           ) 
    6768      !!---------------------------------------------------------------------------------- 
    6869      !!                      ***  ROUTINE  turb_ecmwf  *** 
    … …  
    112113      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    113114      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
     115      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114116      ! 
    115117      INTEGER :: j_itt 
    … …  
    266268            dt_zu = t_zu - sst ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 
    267269            dq_zu = q_zu - ssq ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 
     270 
    268271         END IF 
    269272 
    … …  
    271274         ztmp1 = zu + z0 
    272275         ztmp0 = ztmp1*Linv 
    273          func_m = log(ztmp1) - LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0*Linv) 
     276         func_m = log(ztmp1) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 
    274277         func_h = log(ztmp1) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 
    275278 
    … …  
    280283      ztmp1 = log((zu + z0)/z0q) - psi_h_ecmwf((zu + z0)*Linv) + psi_h_ecmwf(z0q*Linv)   ! func_q 
    281284      Ce = vkarmn*vkarmn/(func_m*ztmp1) 
     285 
     286      ztmp1 = zu + z0 
     287      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 
     288      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 
     289      Cen = vkarmn*vkarmn / (log(ztmp1/z0q)*log(ztmp1/z0q)) 
    282290 
    283291      CALL wrk_dealloc( jpi,jpj,   u_star, t_star, q_star, func_m, func_h, dt_zu, dq_zu, Linv ) 

  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_ncar.F90

    r7753 r8585  
    4848   !                              ! NCAR own values for given constants: 
    4949   REAL(wp), PARAMETER ::   rctv0 = 0.608   ! constant to obtain virtual temperature... 
    50     
    5150   !!---------------------------------------------------------------------- 
    5251CONTAINS 
    5352 
    5453   SUBROUTINE turb_ncar( zt, zu, sst, t_zt, ssq, q_zt, U_zu, & 
    55       &                  Cd, Ch, Ce, t_zu, q_zu, U_blk ) 
     54      &                  Cd, Ch, Ce, t_zu, q_zu, U_blk,      & 
     55      &                  Cdn, Chn, Cen                       ) 
     56 
    5657      !!---------------------------------------------------------------------------------- 
    5758      !!                      ***  ROUTINE  turb_ncar  *** 
    … …  
    112113      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg] 
    113114      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind at 10m                          [m/s] 
     115      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients 
    114116      ! 
    115117      INTEGER ::   j_itt 
    … …  
    199201            ztmp0 = MAX( 0.25 , U_blk/(1. + sqrt_Cd_n10/vkarmn*(LOG(zu/10.) - ztmp2)) ) ! U_n10 (ztmp2 == psi_m(zeta_u)) 
    200202            ztmp0 = cd_neutral_10m(ztmp0)                                               ! Cd_n10 
     203            Cdn(:,:) = ztmp0 
    201204            sqrt_Cd_n10 = sqrt(ztmp0) 
    202205 
    203206            stab    = 0.5 + sign(0.5,zeta_u)                           ! update stability 
    204207            Cx_n10  = 1.e-3*sqrt_Cd_n10*(18.*stab + 32.7*(1. - stab))  ! L&Y 2004 eq. (6c-6d)    (Cx_n10 == Ch_n10) 
     208            Chn(:,:) = Cx_n10 
    205209 
    206210            !! Update of transfer coefficients: 
    … …  
    216220 
    217221         Cx_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)  ! L&Y 2004 eq. (6b)    ! Cx_n10 == Ce_n10 
     222         Cen(:,:) = Cx_n10 
    218223         ztmp1 = 1. + Cx_n10*ztmp0 
    219224         Ce  = Cx_n10*ztmp2 / ztmp1  ! L&Y 2004 eq. (10c) 
Note: See TracChangeset for help on using the changeset viewer.