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 6763 for branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC – NEMO

Ignore:
Timestamp:
2016-06-30T17:17:35+02:00 (8 years ago)
Author:
clem
Message:

new parameterization from Lupkes et al. (2012) of ice-atm and ocean-atm drags (dependent on ice morphology) => ln_Cd_L12 in namelist_ref

Location:
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r6399 r6763  
    1616   !!            3.4  !  2011-11  (C. Harris) Fill arrays required by CICE 
    1717   !!            3.7  !  2014-06  (L. Brodeau) simplification and optimization of CORE bulk 
     18   !!            4.0  !  2016-06  (C. Rousset) Add new param of drags with sea-ice (Lupkes at al 2012) 
    1819   !!---------------------------------------------------------------------- 
    1920 
     
    4546   USE lib_fortran     ! to use key_nosignedzero 
    4647#if defined key_lim3 
    47    USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b 
     48   USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b, at_i_b 
    4849   USE limthd_dh       ! for CALL lim_thd_snwblow 
    4950#elif defined key_lim2 
     
    6061   PUBLIC   blk_ice_core_flx     ! routine called in sbc_ice_lim module 
    6162#endif 
    62    PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module 
     63   PUBLIC   turb_core_2z         ! routine called in sbcblk_mfs module 
    6364 
    6465   INTEGER , PARAMETER ::   jpfld   = 9           ! maximum number of files to read 
     
    7677 
    7778   !                                             !!! CORE bulk parameters 
    78    REAL(wp), PARAMETER ::   rhoa =    1.22        ! air density 
    79    REAL(wp), PARAMETER ::   cpa  = 1000.5         ! specific heat of air 
    80    REAL(wp), PARAMETER ::   Lv   =    2.5e6       ! latent heat of vaporization 
    81    REAL(wp), PARAMETER ::   Ls   =    2.839e6     ! latent heat of sublimation 
    82    REAL(wp), PARAMETER ::   Stef =    5.67e-8     ! Stefan Boltzmann constant 
    83    REAL(wp), PARAMETER ::   Cice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice 
    84    REAL(wp), PARAMETER ::   albo =    0.066       ! ocean albedo assumed to be constant 
     79   REAL(wp), PARAMETER ::   rhoa   =    1.22        ! air density 
     80   REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air 
     81   REAL(wp), PARAMETER ::   Lv     =    2.5e6       ! latent heat of vaporization 
     82   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation 
     83   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant 
     84   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice 
     85   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant 
    8586 
    8687   !                                  !!* Namelist namsbc_core : CORE bulk parameters 
     
    9192   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements 
    9293   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements 
    93  
     94   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012) 
     95 
     96   ! 
     97   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
     98    
    9499   !! * Substitutions 
    95100#  include "domzgr_substitute.h90" 
     
    102107CONTAINS 
    103108 
     109   INTEGER FUNCTION sbc_blk_core_alloc() 
     110      !!------------------------------------------------------------------- 
     111      !!             ***  ROUTINE sbc_blk_core_alloc (clem) *** 
     112      !!------------------------------------------------------------------- 
     113      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_core_alloc ) 
     114      ! 
     115      IF( lk_mpp                  )   CALL mpp_sum( sbc_blk_core_alloc ) 
     116      IF( sbc_blk_core_alloc /= 0 )   CALL ctl_warn('sbc_blk_core_alloc: failed to allocate arrays') 
     117   END FUNCTION sbc_blk_core_alloc 
     118 
     119    
    104120   SUBROUTINE sbc_blk_core( kt ) 
    105121      !!--------------------------------------------------------------------- 
     
    151167         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
    152168         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           & 
    153          &                  sn_tdif, rn_zqt,  rn_zu 
     169         &                  sn_tdif, rn_zqt,  rn_zu    , ln_Cd_L12 
    154170      !!--------------------------------------------------------------------- 
    155171      ! 
     
    157173      IF( kt == nit000 ) THEN                   !  First call kt=nit000  ! 
    158174         !                                      ! ====================== ! 
     175         ! 
     176         !                                      ! allocate sbc_blk_core array (clem) 
     177         IF( sbc_blk_core_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_core : unable to allocate standard arrays' ) 
    159178         ! 
    160179         REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters 
     
    319338         &               Cd, Ch, Ce, zt_zu, zq_zu ) 
    320339     
     340      ! Make ocean-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 
     341#if defined key_lim3 
     342      IF( ln_Cd_L12 ) THEN 
     343 
     344         Cd_oce(:,:) = Cd(:,:)    ! record value of pure ocean-atm. drag 
     345 
     346         CALL Cdn10_Lupkes2012( Cd )  ! calculate new drag from Lupkes(2012) equations 
     347 
     348      ENDIF 
     349#endif 
     350       
    321351      ! ... tau module, i and j component 
    322352      DO jj = 1, jpj 
     
    437467      !!--------------------------------------------------------------------- 
    438468      INTEGER  ::   ji, jj    ! dummy loop indices 
    439       REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2 
    440469      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point 
    441470      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point 
     471      REAL(wp), DIMENSION(:,:), POINTER ::   Cd               ! transfer coefficient for momentum      (tau) 
    442472      !!--------------------------------------------------------------------- 
    443473      ! 
    444474      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_tau') 
    445475      ! 
    446       ! local scalars ( place there for vector optimisation purposes) 
    447       zcoef_wnorm  = rhoa * Cice 
    448       zcoef_wnorm2 = rhoa * Cice * 0.5 
     476      CALL wrk_alloc( jpi,jpj, Cd ) 
     477 
     478      Cd(:,:) = Cd_ice 
     479       
     480      ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 
     481#if defined key_lim3 
     482      IF( ln_Cd_L12 ) THEN 
     483 
     484         CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 
     485 
     486      ENDIF 
     487#endif 
    449488 
    450489!!gm brutal.... 
     
    467506               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    468507                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    469                zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
     508               zwnorm_f = rhoa * Cd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    470509               ! ... ice stress at I-point 
    471510               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
     
    476515               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
    477516                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
    478                wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     517               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    479518            END DO 
    480519         END DO 
     
    493532         DO jj = 2, jpjm1 
    494533            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    495                utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     534               utau_ice(ji,jj) = rhoa * Cd(ji,jj) * 0.5_wp * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
    496535                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    497                vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     536               vtau_ice(ji,jj) = rhoa * Cd(ji,jj) * 0.5_wp * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
    498537                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    499538            END DO 
     
    510549      ENDIF 
    511550 
     551      CALL wrk_dealloc( jpi,jpj, Cd ) 
     552 
    512553      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau') 
    513554       
     
    548589      ! local scalars ( place there for vector optimisation purposes) 
    549590      zcoef_dqlw   = 4.0 * 0.95 * Stef 
    550       zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
    551       zcoef_dqsb   = rhoa * cpa * Cice 
     591      zcoef_dqla   = -Ls * Cd_ice * 11637800. * (-5897.8) 
     592      zcoef_dqsb   = rhoa * cpa * Cd_ice 
    552593 
    553594      zztmp = 1. / ( 1. - albo ) 
     
    575616               ! ... turbulent heat fluxes 
    576617               ! Sensible Heat 
    577                z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
     618               z_qsb(ji,jj,jl) = rhoa * cpa * Cd_ice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    578619               ! Latent Heat 
    579                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                            
     620               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cd_ice * wndm_ice(ji,jj)   &                            
    580621                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    581622              ! Latent heat sensitivity for ice (Dqla/Dt) 
     
    820861         ! 
    821862      END DO 
    822  
     863       
    823864      CALL wrk_dealloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 
    824865      CALL wrk_dealloc( jpi,jpj, zeta_u, stab ) 
     
    903944   END FUNCTION psi_h 
    904945 
     946 
     947   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
     948      !!---------------------------------------------------------------------- 
     949      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
     950      !! 
     951      !! ** Purpose :    Recompute the ice-atm and ocean-atm drags at 10m height to make 
     952      !!                 them dependent on edges at leads, melt ponds and flows. 
     953      !!                 After some approximations, this can be resumed to a dependency 
     954      !!                 on ice concentration. 
     955      !!                 
     956      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50) 
     957      !!                 with the highest level of approximation: level4, eq.(59) 
     958      !!                 The drag can be re-written as follows: 
     959      !! 
     960      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu 
     961      !! 
     962      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59) 
     963      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59) 
     964      !!                    A is the concentration of ice minus melt ponds (if any) 
     965      !! 
     966      !!                 This new drag has a parabolic shape (as a function of A) starting at 
     967      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5  
     968      !!                 and going down to Cdi(say 1.4e-3) for A=1 
     969      !! 
     970      !!                 It is theoretically applicable to all ice conditions (not only MIZ) 
     971      !!                 => see Lupkes et al (2013) 
     972      !! 
     973      !! ** References : Lupkes et al. JGR 2012 (theory) 
     974      !!                 Lupkes et al. GRL 2013 (application to GCM) 
     975      !! 
     976      !!---------------------------------------------------------------------- 
     977      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
     978      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp 
     979      REAL(wp), PARAMETER ::   znu   = 1._wp 
     980      REAL(wp), PARAMETER ::   zmu   = 1._wp 
     981      REAL(wp), PARAMETER ::   zbeta = 1._wp 
     982      REAL(wp)            ::   zcoef 
     983      !!---------------------------------------------------------------------- 
     984      zcoef = znu + 1._wp / ( 10._wp * zbeta ) 
     985 
     986      Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag 
     987         &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag 
     988         &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology 
     989       
     990   END SUBROUTINE Cdn10_Lupkes2012 
     991       
    905992   !!====================================================================== 
    906993END MODULE sbcblk_core 
  • branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r6746 r6763  
    617617      u_ice_b(:,:)     = u_ice(:,:) 
    618618      v_ice_b(:,:)     = v_ice(:,:) 
     619      at_i_b (:,:)     = SUM( a_i_b(:,:,:), dim=3 ) 
    619620       
    620621   END SUBROUTINE sbc_lim_bef 
Note: See TracChangeset for help on using the changeset viewer.