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

Ignore:
Timestamp:
2017-11-24T17:56:51+01:00 (6 years ago)
Author:
gm
Message:

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

File:
1 edited

Legend:

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

    r8637 r8813  
    1515   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk 
    1616   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore 
    17    !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 
     17   !!                 !                        ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 
    1818   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine 
     19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle)  
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    2324   !!   sbc_blk       : bulk formulation as ocean surface boundary condition 
    2425   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean 
    25    !!   blk_ice       : computes momentum, heat and freshwater fluxes over sea ice 
    2626   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP 
    2727   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air) 
    2828   !!   q_sat         : saturation humidity as a function of SLP and temperature 
    2929   !!   L_vap         : latent heat of vaporization of water as a function of temperature 
     30   !!             sea-ice case only :  
     31   !!   blk_ice_tau   : provide the air-ice stress 
     32   !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface 
     33   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler) 
     34   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
     35   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
    3036   !!---------------------------------------------------------------------- 
    3137   USE oce            ! ocean dynamics and tracers 
     
    4248   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 
    4349   USE icethd_dh      ! for CALL ice_thd_snwblow 
     50   USE icethd_zdf,ONLY:   rn_cnd_s  ! for blk_ice_qcn 
    4451#endif 
    4552   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
     
    6370   PUBLIC   blk_ice_tau   ! routine called in icestp module 
    6471   PUBLIC   blk_ice_flx   ! routine called in icestp module 
    65 #endif 
     72   PUBLIC   blk_ice_qcn   ! routine called in icestp module 
     73#endif  
    6674 
    6775!!Lolo: should ultimately be moved in the module with all physical constants ? 
     
    111119   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 
    112120   ! 
    113 !!cr   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
    114121   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau) 
    115122   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens) 
     
    139146      !!             ***  ROUTINE sbc_blk_alloc *** 
    140147      !!------------------------------------------------------------------- 
    141 !!cr      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 
    142148      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 
    143149         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 
     
    147153   END FUNCTION sbc_blk_alloc 
    148154 
     155 
    149156   SUBROUTINE sbc_blk_init 
    150157      !!--------------------------------------------------------------------- 
     
    154161      !! 
    155162      !! ** Method  :  
    156       !! 
    157       !!      C A U T I O N : never mask the surface stress fields 
    158       !!                      the stress is assumed to be in the (i,j) mesh referential 
    159       !! 
    160       !! ** Action  :    
    161163      !! 
    162164      !!---------------------------------------------------------------------- 
     
    285287      !! 
    286288      !! ** Purpose :   provide at each time step the surface ocean fluxes 
    287       !!      (momentum, heat, freshwater and runoff) 
     289      !!              (momentum, heat, freshwater and runoff) 
    288290      !! 
    289291      !! ** Method  : (1) READ each fluxes in NetCDF files: 
     
    566568   END SUBROUTINE blk_oce 
    567569 
     570 
     571 
     572   FUNCTION rho_air( ptak, pqa, pslp ) 
     573      !!------------------------------------------------------------------------------- 
     574      !!                           ***  FUNCTION rho_air  *** 
     575      !! 
     576      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
     577      !! 
     578      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
     579      !!------------------------------------------------------------------------------- 
     580      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
     581      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
     582      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
     583      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
     584      !!------------------------------------------------------------------------------- 
     585      ! 
     586      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
     587      ! 
     588   END FUNCTION rho_air 
     589 
     590 
     591   FUNCTION cp_air( pqa ) 
     592      !!------------------------------------------------------------------------------- 
     593      !!                           ***  FUNCTION cp_air  *** 
     594      !! 
     595      !! ** Purpose : Compute specific heat (Cp) of moist air 
     596      !! 
     597      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     598      !!------------------------------------------------------------------------------- 
     599      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
     600      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
     601      !!------------------------------------------------------------------------------- 
     602      ! 
     603      Cp_air = Cp_dry + Cp_vap * pqa 
     604      ! 
     605   END FUNCTION cp_air 
     606 
     607 
     608   FUNCTION q_sat( ptak, pslp ) 
     609      !!---------------------------------------------------------------------------------- 
     610      !!                           ***  FUNCTION q_sat  *** 
     611      !! 
     612      !! ** Purpose : Specific humidity at saturation in [kg/kg] 
     613      !!              Based on accurate estimate of "e_sat" 
     614      !!              aka saturation water vapor (Goff, 1957) 
     615      !! 
     616      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     617      !!---------------------------------------------------------------------------------- 
     618      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K] 
     619      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa] 
     620      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg] 
     621      ! 
     622      INTEGER  ::   ji, jj         ! dummy loop indices 
     623      REAL(wp) ::   ze_sat, ztmp   ! local scalar 
     624      !!---------------------------------------------------------------------------------- 
     625      ! 
     626      DO jj = 1, jpj 
     627         DO ji = 1, jpi 
     628            ! 
     629            ztmp = rt0 / ptak(ji,jj) 
     630            ! 
     631            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
     632            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        & 
     633               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  & 
     634               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
     635               ! 
     636            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa] 
     637            ! 
     638         END DO 
     639      END DO 
     640      ! 
     641   END FUNCTION q_sat 
     642 
     643 
     644   FUNCTION gamma_moist( ptak, pqa ) 
     645      !!---------------------------------------------------------------------------------- 
     646      !!                           ***  FUNCTION gamma_moist  *** 
     647      !! 
     648      !! ** Purpose : Compute the moist adiabatic lapse-rate. 
     649      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
     650      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
     651      !! 
     652      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     653      !!---------------------------------------------------------------------------------- 
     654      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
     655      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
     656      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
     657      ! 
     658      INTEGER  ::   ji, jj         ! dummy loop indices 
     659      REAL(wp) :: zrv, ziRT        ! local scalar 
     660      !!---------------------------------------------------------------------------------- 
     661      ! 
     662      DO jj = 1, jpj 
     663         DO ji = 1, jpi 
     664            zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
     665            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
     666            gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
     667         END DO 
     668      END DO 
     669      ! 
     670   END FUNCTION gamma_moist 
     671 
     672 
     673   FUNCTION L_vap( psst ) 
     674      !!--------------------------------------------------------------------------------- 
     675      !!                           ***  FUNCTION L_vap  *** 
     676      !! 
     677      !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
     678      !! 
     679      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     680      !!---------------------------------------------------------------------------------- 
     681      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
     682      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
     683      !!---------------------------------------------------------------------------------- 
     684      ! 
     685      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
     686      ! 
     687   END FUNCTION L_vap 
     688 
    568689#if defined key_lim3 
     690   !!---------------------------------------------------------------------- 
     691   !!   'key_lim3'                                       ESIM sea-ice model 
     692   !!---------------------------------------------------------------------- 
     693   !!   blk_ice_tau : provide the air-ice stress 
     694   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface 
     695   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler) 
     696   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
     697   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
     698   !!---------------------------------------------------------------------- 
    569699 
    570700   SUBROUTINE blk_ice_tau 
     
    688818 
    689819 
    690    SUBROUTINE blk_ice_flx( ptsu, palb ) 
     820   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 
    691821      !!--------------------------------------------------------------------- 
    692822      !!                     ***  ROUTINE blk_ice_flx  *** 
    693823      !! 
    694       !! ** Purpose :   provide the surface boundary condition over sea-ice 
     824      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface 
    695825      !! 
    696826      !! ** Method  :   compute heat and freshwater exchanged 
     
    701831      !!--------------------------------------------------------------------- 
    702832      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature 
     833      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness 
     834      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness 
    703835      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies) 
    704836      !! 
     
    707839      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    708840      REAL(wp) ::   zztmp, z1_lsub           !   -      - 
     841      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice 
     842      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables 
    709843      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice 
    710844      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice 
     
    717851      IF( ln_timing )   CALL timing_start('blk_ice_flx') 
    718852      ! 
    719       ! 
    720       ! local scalars 
    721       zcoef_dqlw   = 4.0 * 0.95 * Stef 
    722       zcoef_dqla   = -Ls * 11637800. * (-5897.8) 
     853      zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars 
     854      zcoef_dqla = -Ls * 11637800. * (-5897.8) 
    723855      ! 
    724856      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     
    815947      END DO 
    816948 
    817       !-------------------------------------------------------------------- 
    818       ! FRACTIONs of net shortwave radiation which is not absorbed in the 
    819       ! thin surface layer and penetrates inside the ice cover 
    820       ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    821       ! 
    822       fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    823       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     949      ! --- absorbed and transmitted shortwave radiation (W/m2) --- ! 
     950      ! 
     951      ! former coding was 
     952      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     953      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     954 
     955      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- ! 
     956      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value 
     957      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1 
     958 
     959      qsr_ice_tr(:,:,:) = 0._wp 
     960 
     961      DO jl = 1, jpl 
     962         DO jj = 1, jpj 
     963            DO ji = 1, jpi 
     964               ! 
     965               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1 
     966               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 
     967               ! 
     968               IF ( phs(ji,jj,jl) <= 0._wp ) THEN   ;    zfrqsr_tr_i  = zfr1 + zfac * zfr2 
     969               ELSE                                 ;    zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque 
     970               ENDIF 
     971               ! 
     972               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation  
     973               ! 
     974            END DO 
     975         END DO 
     976      END DO 
    824977      ! 
    825978      IF(ln_ctl) THEN 
     
    836989   END SUBROUTINE blk_ice_flx 
    837990    
    838 #endif 
    839  
    840    FUNCTION rho_air( ptak, pqa, pslp ) 
    841       !!------------------------------------------------------------------------------- 
    842       !!                           ***  FUNCTION rho_air  *** 
    843       !! 
    844       !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    845       !! 
    846       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
    847       !!------------------------------------------------------------------------------- 
    848       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
    849       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
    850       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
    851       REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
    852       !!------------------------------------------------------------------------------- 
    853       ! 
    854       rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
    855       ! 
    856    END FUNCTION rho_air 
    857  
    858  
    859    FUNCTION cp_air( pqa ) 
    860       !!------------------------------------------------------------------------------- 
    861       !!                           ***  FUNCTION cp_air  *** 
    862       !! 
    863       !! ** Purpose : Compute specific heat (Cp) of moist air 
    864       !! 
    865       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    866       !!------------------------------------------------------------------------------- 
    867       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    868       REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
    869       !!------------------------------------------------------------------------------- 
    870       ! 
    871       Cp_air = Cp_dry + Cp_vap * pqa 
    872       ! 
    873    END FUNCTION cp_air 
    874  
    875  
    876    FUNCTION q_sat( ptak, pslp ) 
    877       !!---------------------------------------------------------------------------------- 
    878       !!                           ***  FUNCTION q_sat  *** 
    879       !! 
    880       !! ** Purpose : Specific humidity at saturation in [kg/kg] 
    881       !!              Based on accurate estimate of "e_sat" 
    882       !!              aka saturation water vapor (Goff, 1957) 
    883       !! 
    884       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    885       !!---------------------------------------------------------------------------------- 
    886       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K] 
    887       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa] 
    888       REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg] 
    889       ! 
    890       INTEGER  ::   ji, jj         ! dummy loop indices 
    891       REAL(wp) ::   ze_sat, ztmp   ! local scalar 
    892       !!---------------------------------------------------------------------------------- 
    893       ! 
    894       DO jj = 1, jpj 
    895          DO ji = 1, jpi 
    896             ! 
    897             ztmp = rt0 / ptak(ji,jj) 
    898             ! 
    899             ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
    900             ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        & 
    901                &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  & 
    902                &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
     991 
     992   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi ) 
     993      !!--------------------------------------------------------------------- 
     994      !!                     ***  ROUTINE blk_ice_qcn  *** 
     995      !! 
     996      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux 
     997      !!                to force sea ice / snow thermodynamics 
     998      !!                in the case JULES coupler is emulated 
     999      !!                 
     1000      !! ** Method  :   compute surface energy balance assuming neglecting heat storage 
     1001      !!                following the 0-layer Semtner (1976) approach 
     1002      !! 
     1003      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K) 
     1004      !!              - qcn_ice : surface inner conduction flux (W/m2) 
     1005      !! 
     1006      !!--------------------------------------------------------------------- 
     1007      INTEGER                   , INTENT(in   ) ::   k_monocat   ! single-category option 
     1008      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu        ! sea ice / snow surface temperature 
     1009      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb         ! sea ice base temperature 
     1010      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs         ! snow thickness 
     1011      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi         ! sea ice thickness 
     1012      ! 
     1013      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations 
     1014      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction 
     1015      ! 
     1016      INTEGER  ::   ji, jj, jl           ! dummy loop indices 
     1017      INTEGER  ::   iter                 ! local integer 
     1018      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars 
     1019      REAL(wp) ::   zkeff_h, ztsu        ! 
     1020      REAL(wp) ::   zqc, zqnet           ! 
     1021      REAL(wp) ::   zhe, zqa0            ! 
     1022      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor 
     1023      !!--------------------------------------------------------------------- 
     1024 
     1025      IF( nn_timing == 1 )  CALL timing_start('blk_ice_qcn') 
     1026       
     1027      ! -------------------------------------! 
     1028      !      I   Enhanced conduction factor  ! 
     1029      ! -------------------------------------! 
     1030      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3) 
     1031      ! Fichefet and Morales Maqueda, JGR 1997 
     1032      ! 
     1033      zgfac(:,:,:) = 1._wp 
     1034       
     1035      SELECT CASE ( k_monocat )  
     1036      ! 
     1037      CASE ( 1 , 3 ) 
     1038         ! 
     1039         zfac    =   1._wp /  ( rn_cnd_s + rcdic ) 
     1040         zfac2   =   EXP(1._wp) * 0.5_wp * zepsilon 
     1041         zfac3   =   2._wp / zepsilon 
     1042         !    
     1043         DO jl = 1, jpl                 
     1044            DO jj = 1 , jpj 
     1045               DO ji = 1, jpi 
     1046                  !                                ! Effective thickness 
     1047                  zhe               =   ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac 
     1048                  !                                ! Enhanced conduction factor 
     1049                  IF( zhe >=  zfac2 )  & 
     1050                     zgfac(ji,jj,jl)   =   MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) ) 
     1051               END DO 
     1052            END DO 
     1053         END DO 
     1054         !       
     1055      END SELECT 
     1056       
     1057      ! -------------------------------------------------------------! 
     1058      !      II   Surface temperature and conduction flux            ! 
     1059      ! -------------------------------------------------------------! 
     1060      ! 
     1061      zfac   =   rcdic * rn_cnd_s  
     1062      !                             ! ========================== ! 
     1063      DO jl = 1, jpl                !  Loop over ice categories  ! 
     1064         !                          ! ========================== ! 
     1065         DO jj = 1 , jpj 
     1066            DO ji = 1, jpi 
     1067               !                    ! Effective conductivity of the snow-ice system divided by thickness 
     1068               zkeff_h =   zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) ) 
     1069               !                    ! Store initial surface temperature 
     1070               ztsu    =   ptsu(ji,jj,jl) 
     1071               !                    ! Net initial atmospheric heat flux 
     1072               zqa0    =   qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl) 
    9031073               ! 
    904             q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa] 
    905             ! 
    906          END DO 
    907       END DO 
    908       ! 
    909    END FUNCTION q_sat 
    910  
    911  
    912    FUNCTION gamma_moist( ptak, pqa ) 
    913       !!---------------------------------------------------------------------------------- 
    914       !!                           ***  FUNCTION gamma_moist  *** 
    915       !! 
    916       !! ** Purpose : Compute the moist adiabatic lapse-rate. 
    917       !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
    918       !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
    919       !! 
    920       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    921       !!---------------------------------------------------------------------------------- 
    922       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
    923       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
    924       REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
    925       ! 
    926       INTEGER  ::   ji, jj         ! dummy loop indices 
    927       REAL(wp) :: zrv, ziRT        ! local scalar 
    928       !!---------------------------------------------------------------------------------- 
    929       ! 
    930       DO jj = 1, jpj 
    931          DO ji = 1, jpi 
    932             zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
    933             ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
    934             gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
    935          END DO 
    936       END DO 
    937       ! 
    938    END FUNCTION gamma_moist 
    939  
    940  
    941    FUNCTION L_vap( psst ) 
    942       !!--------------------------------------------------------------------------------- 
    943       !!                           ***  FUNCTION L_vap  *** 
    944       !! 
    945       !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
    946       !! 
    947       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    948       !!---------------------------------------------------------------------------------- 
    949       REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
    950       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    951       !!---------------------------------------------------------------------------------- 
    952       ! 
    953       L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
    954       ! 
    955    END FUNCTION L_vap 
    956  
    957 #if defined key_lim3 
     1074               DO iter = 1, nit     ! --- Iteration loop 
     1075                   !                      ! Conduction heat flux through snow-ice system (>0 downwards) 
     1076                   zqc   =   zkeff_h * ( ztsu - ptb(ji,jj) ) 
     1077                   !                      ! Surface energy budget 
     1078                   zqnet =   zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc 
     1079                   !                      ! Temperature update 
     1080                   ztsu  =   ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) 
     1081               END DO 
     1082               ! 
     1083               ptsu   (ji,jj,jl) = MIN( rt0, ztsu ) 
     1084               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 
     1085               ! 
     1086            END DO 
     1087         END DO 
     1088         ! 
     1089      END DO  
     1090      !       
     1091      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_qcn') 
     1092      ! 
     1093   END SUBROUTINE blk_ice_qcn 
     1094    
    9581095 
    9591096   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
     
    9611098      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
    9621099      !! 
    963       !! ** Purpose :    Recompute the ice-atm drag at 10m height to make 
    964       !!                 it dependent on edges at leads, melt ponds and flows. 
     1100      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m  
     1101      !!                 to make it dependent on edges at leads, melt ponds and flows. 
    9651102      !!                 After some approximations, this can be resumed to a dependency 
    9661103      !!                 on ice concentration. 
     
    10121149      !!                      ***  ROUTINE  Cdn10_Lupkes2015  *** 
    10131150      !! 
    1014       !! ** pUrpose :    1lternative turbulent transfert coefficients formulation 
     1151      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation 
    10151152      !!                 between sea-ice and atmosphere with distinct momentum  
    10161153      !!                 and heat coefficients depending on sea-ice concentration  
Note: See TracChangeset for help on using the changeset viewer.