Changeset 758


Ignore:
Timestamp:
2007-12-07T15:44:32+01:00 (14 years ago)
Author:
ctlod
Message:

Use the previous version TURB_CORE_1Z from the reference which keep same results, see ticket:#38

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/dev_001_SBC/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r710 r758  
    174174      !                                               ! input fields at the current time-step 
    175175 
    176 !!gm  IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    177  
    178       CALL blk_oce_core( sst_m, ssu_m, ssv_m )        ! set the ocean surface fluxes 
    179  
    180 !!gm  ENDIF 
     176      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
     177 
     178          CALL blk_oce_core( sst_m, ssu_m, ssv_m )        ! set the ocean surface fluxes 
     179 
     180      ENDIF 
    181181      !                                               ! using CORE bulk formulea 
    182182   END SUBROUTINE sbc_blk_core 
     
    271271!        &                 Cd(:,:), Ch(:,:), Ce(:,:)                              ) 
    272272!!gm end 
    273       CALL TURB_CORE( 10., zst   , sf(jp_tair)%fnow,   & 
    274          &                 zqsatw, sf(jp_humi)%fnow, zwind_speed_t,   & 
    275          &                 Cd, Ch, Ce                              ) 
     273      !! If air temp. and spec. hum. are given at same height than wind (10m) : 
     274      CALL TURB_CORE_1Z(10., zst   , sf(jp_tair)%fnow,                  & 
     275         &                   zqsatw, sf(jp_humi)%fnow,  zwind_speed_t,  & 
     276         &                   Cd, Ch, Ce ) 
    276277    
    277278      ! ...  umasked Momentum : utau, vtau at U- and V_points, resp. 
     
    465466   END SUBROUTINE blk_ice_core 
    466467   
     468   SUBROUTINE TURB_CORE_1Z(zu, sst, T_a, q_sat, q_a,   & 
     469      &                        dU, Cd, Ch, Ce   ) 
     470      !!---------------------------------------------------------------------- 
     471      !!                      ***  ROUTINE  turb_core  *** 
     472      !! 
     473      !! ** Purpose :   Computes turbulent transfert coefficients of surface 
     474      !!                fluxes according to Large & Yeager (2004) 
     475      !! 
     476      !! ** Method  :   I N E R T I A L   D I S S I P A T I O N   M E T H O D 
     477      !!      Momentum, Latent and sensible heat exchange coefficients 
     478      !!      Caution: this procedure should only be used in cases when air 
     479      !!      temperature (T_air), air specific humidity (q_air) and wind (dU) 
     480      !!      are provided at the same height 'zzu'! 
     481      !! 
     482      !! References : 
     483      !!      Large & Yeager, 2004 : ??? 
     484      !! History : 
     485      !!        !  XX-XX  (???     )  Original code 
     486      !!   9.0  !  05-08  (L. Brodeau) Rewriting and optimization 
     487      !!---------------------------------------------------------------------- 
     488      !! * Arguments 
     489 
     490      REAL(wp), INTENT(in) :: zu                 ! altitude of wind measurement       [m] 
     491      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::  & 
     492         sst,       &       ! sea surface temperature         [Kelvin] 
     493         T_a,       &       ! potential air temperature       [Kelvin] 
     494         q_sat,     &       ! sea surface specific humidity   [kg/kg] 
     495         q_a,       &       ! specific air humidity           [kg/kg] 
     496         dU                 ! wind module |U(zu)-U(0)|        [m/s] 
     497      REAL(wp), intent(out), DIMENSION(jpi,jpj)  :: & 
     498         Cd,    &                ! transfert coefficient for momentum       (tau) 
     499         Ch,    &                ! transfert coefficient for temperature (Q_sens) 
     500         Ce                      ! transfert coefficient for evaporation  (Q_lat) 
     501 
     502      !! * Local declarations 
     503      REAL(wp), DIMENSION(jpi,jpj)  ::   & 
     504         dU10,        &       ! dU                                   [m/s] 
     505         dT,          &       ! air/sea temperature differeence      [K] 
     506         dq,          &       ! air/sea humidity difference          [K] 
     507         Cd_n10,      &       ! 10m neutral drag coefficient 
     508         Ce_n10,      &       ! 10m neutral latent coefficient 
     509         Ch_n10,      &       ! 10m neutral sensible coefficient 
     510         sqrt_Cd_n10, &       ! root square of Cd_n10 
     511         sqrt_Cd,     &       ! root square of Cd 
     512         T_vpot,      &       ! virtual potential temperature        [K] 
     513         T_star,      &       ! turbulent scale of tem. fluct. 
     514         q_star,      &       ! turbulent humidity of temp. fluct. 
     515         U_star,      &       ! turb. scale of velocity fluct. 
     516         L,           &       ! Monin-Obukov length                  [m] 
     517         zeta,        &       ! stability parameter at height zu 
     518         U_n10,       &       ! neutral wind velocity at 10m         [m]    
     519         xlogt, xct, zpsi_h, zpsi_m 
     520      !! 
     521      INTEGER :: j_itt 
     522      INTEGER, PARAMETER :: nb_itt = 3 
     523      INTEGER, DIMENSION(jpi,jpj)  ::   & 
     524         stab         ! 1st guess stability test integer 
     525 
     526      REAL(wp), PARAMETER ::                        & 
     527         grav   = 9.8,          &  ! gravity                        
     528         kappa  = 0.4              ! von Karman s constant 
     529 
     530      !! * Start 
     531      !! Air/sea differences 
     532      dU10 = max(0.5, dU)     ! we don't want to fall under 0.5 m/s 
     533      dT = T_a - sst          ! assuming that T_a is allready the potential temp. at zzu 
     534      dq = q_a - q_sat 
     535      !!     
     536      !! Virtual potential temperature 
     537      T_vpot = T_a*(1. + 0.608*q_a) 
     538      !! 
     539      !! Neutral Drag Coefficient 
     540      stab    = 0.5 + sign(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0  
     541      Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a) 
     542      sqrt_Cd_n10 = sqrt(Cd_n10) 
     543      Ce_n10  = 1E-3 * ( 34.6 * sqrt_Cd_n10 )               !   L & Y eq. (6b) 
     544      Ch_n10  = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1-stab)) !   L & Y eq. (6c), (6d) 
     545      !! 
     546      !! Initializing transfert coefficients with their first guess neutral equivalents : 
     547      Cd = Cd_n10 ;  Ce = Ce_n10 ;  Ch = Ch_n10 ;  sqrt_Cd = sqrt(Cd) 
     548 
     549      !! * Now starting iteration loop 
     550      DO j_itt=1, nb_itt 
     551         !! Turbulent scales : 
     552         U_star  = sqrt_Cd*dU10                !   L & Y eq. (7a) 
     553         T_star  = Ch/sqrt_Cd*dT               !   L & Y eq. (7b) 
     554         q_star  = Ce/sqrt_Cd*dq               !   L & Y eq. (7c) 
     555 
     556         !! Estimate the Monin-Obukov length : 
     557         L  = (U_star**2)/( kappa*grav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) ) 
     558 
     559         !! Stability parameters : 
     560         zeta = zu/L ;  zeta   = sign( min(abs(zeta),10.0), zeta ) 
     561         zpsi_h = psi_h(zeta) 
     562         zpsi_m = psi_m(zeta) 
     563 
     564         !! Shifting the wind speed to 10m and neutral stability : 
     565         U_n10 = dU10*1./(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) !  L & Y eq. (9a) 
     566 
     567         !! Updating the neutral 10m transfer coefficients : 
     568         Cd_n10  = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a) 
     569         sqrt_Cd_n10 = sqrt(Cd_n10) 
     570         Ce_n10  = 1E-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b) 
     571         stab    = 0.5 + sign(0.5,zeta) 
     572         Ch_n10  = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab))           !  L & Y eq. (6c), (6d) 
     573 
     574         !! Shifting the neutral  10m transfer coefficients to ( zu , zeta ) : 
     575         !! 
     576         xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10) - zpsi_m) 
     577         Cd  = Cd_n10/(xct*xct) ;  sqrt_Cd = sqrt(Cd) 
     578         !! 
     579         xlogt = log(zu/10.) - zpsi_h 
     580         !! 
     581         xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10 
     582         Ch  = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct 
     583         !! 
     584         xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10 
     585         Ce  = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct 
     586         !! 
     587      END DO 
     588      !! 
     589    END SUBROUTINE TURB_CORE_1Z 
     590 
     591    FUNCTION psi_m(zta)   !! Psis, L & Y eq. (8c), (8d), (8e) 
     592      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 
     593 
     594      REAL(wp), PARAMETER :: pi = 3.141592653589793_wp 
     595      REAL(wp), DIMENSION(jpi,jpj)             :: psi_m 
     596      REAL(wp), DIMENSION(jpi,jpj)             :: X2, X, stabit 
     597      X2 = sqrt(abs(1. - 16.*zta))  ;  X2 = max(X2 , 1.0) ;  X  = sqrt(X2) 
     598      stabit    = 0.5 + sign(0.5,zta) 
     599      psi_m = -5.*zta*stabit  &                                                  ! Stable 
     600           & + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2)  ! Unstable  
     601    END FUNCTION psi_m 
     602 
     603    FUNCTION psi_h(zta)    !! Psis, L & Y eq. (8c), (8d), (8e) 
     604      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 
     605 
     606      REAL(wp), DIMENSION(jpi,jpj)             :: psi_h 
     607      REAL(wp), DIMENSION(jpi,jpj)             :: X2, X, stabit 
     608      X2 = sqrt(abs(1. - 16.*zta))  ;  X2 = max(X2 , 1.) ;  X  = sqrt(X2) 
     609      stabit    = 0.5 + sign(0.5,zta) 
     610      psi_h = -5.*zta*stabit  &                                       ! Stable 
     611           & + (1. - stabit)*(2.*log( (1. + X2)/2. ))                 ! Unstable 
     612    END FUNCTION psi_h 
    467613   
    468    SUBROUTINE turb_core( zzu, T_0, T_a, q_sat, q_a, dU, C_d, C_h, C_e ) 
    469       !!--------------------------------------------------------------------- 
    470       !!       Computes turbulent transfert coefficients of surface fluxes 
    471       !!             according to Large & Yeager (2004) 
    472       !! 
    473       !!    I N E R T I A L   D I S S I P A T I O N   M E T H O D 
    474       !! 
    475       !!   Momentum, Latent and sensible heat exchange coefficients 
    476       !! 
    477       !!   Caution: this procedure should only be used in cases when air temperature (T_air),  
    478       !!   air specific humidity (q_air) and wind (dU) are provided at the same height 'zzu'! 
    479       !! 
    480       !!   Laurent Brodeau, LEGI, Grenoble 
    481       !!   brodeau@hmg.inpg.fr 
    482       !!--------------------------------------------------------------------- 
    483       REAL(wp), INTENT(in   )                     ::   zzu     ! altitude of wind measurement      [m] 
    484       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   T_0     ! sea surface temperature           [Kelvin] 
    485       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   T_a     ! potential air temperature at zzu  [Kelvin] 
    486       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_sat   ! sea surface specific humidity     [kg/kg] 
    487       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_a     ! specific air humidity at zzu      [kg/kg] 
    488       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   dU      ! wind module |U(zzu)-U(0)|         [m/s] 
    489       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   C_d     ! transfer coefficient for momentum      (tau) 
    490       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   C_h     ! transfer coefficient for sensible heat (Q_sens) 
    491       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   C_e     ! tansfert coefficient for evaporation   (Q_lat) 
    492  
    493       INTEGER             ::   jk                   ! dummy loop indices 
    494       INTEGER , PARAMETER ::   nit = 3              ! number of iterations 
    495       INTEGER , DIMENSION(jpi,jpj) ::   stab        ! stability test integer 
    496       INTEGER , DIMENSION(jpi,jpj) ::   stabit      ! stability within iterative loop 
    497       REAL(wp), DIMENSION(jpi,jpj) ::   & 
    498          dU10,        &  ! dU                                   [m/s] 
    499          dT,          &  ! air/sea temperature differeence      [K] 
    500          dq,          &  ! air/sea humidity difference          [K] 
    501          Cd_n10,      &  ! 10m neutral drag coefficient 
    502          Ce_n10,      &  ! 10m neutral latent coefficient 
    503          Ch_n10,      &  ! 10m neutral sensible coefficient 
    504          Cd,          &  ! drag coefficient 
    505          Ce,          &  ! latent coefficient 
    506          Ch,          &  ! sensible coefficient 
    507          sqrt_Cd_n10, &  ! root square of Cd_n10 
    508          sqrt_Cd,     &  ! root square of Cd 
    509          T_vpot,      &  ! virtual potential temperature        [K] 
    510          T_star,      &  ! turbulent scale of tem. fluct. 
    511          q_star,      &  ! turbulent humidity of temp. fluct. 
    512          U_star,      &  ! turb. scale of velocity fluct. 
    513          L,           &  ! Monin-Obukov length                  [m] 
    514          zeta,        &  ! stability parameter at height zzu 
    515          X2, X,       &    
    516          psi_m,       & 
    517          psi_h,       & 
    518          U_n10,       &  ! neutral wind velocity at 10m        [m]    
    519          xlogt 
    520       !!--------------------------------------------------------------------- 
    521  
    522       !! I. Preliminary stuffs 
    523       !! -------------------- 
    524  
    525       ! ... Air/sea differences 
    526       dU10 = MAX( 0.5, dU )     ! we do not want to fall under 0.5 m/s 
    527       dT   = T_a - T_0          ! assuming that T_a is allready the potential temp. at zzu 
    528       dq   = q_a - q_sat 
    529  
    530       ! ... Virtual potential temperature 
    531       T_vpot = T_a * ( 1. + 0.608 * q_a ) 
    532  
    533       ! ... Computing Neutral Drag Coefficient 
    534 !CDIR NOVERRCHK 
    535       Cd_n10      = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )          !  \\ L & Y eq. (6a) 
    536       sqrt_Cd_n10 = SQRT( Cd_n10 ) 
    537  
    538       Ce_n10  = 1E-3 * ( 34.6 * sqrt_Cd_n10 )                         !  \\ L & Y eq. (6b) 
    539  
    540       ! ... First guess of stabilitty : 
    541       stab    = 0.5 + SIGN( 0.5, dT )     ! stable : stab = 1 ; unstable : stab = 0  
    542       Ch_n10  = 1E-3 * sqrt_Cd_n10 * ( 18*stab + 32.7*(1-stab) )      !  \\ L & Y eq. (6c), (6d) 
    543  
    544       ! ... Initializing transfert coefficients with their first guess neutral equivalents : 
    545       Cd = Cd_n10   ;   Ce = Ce_n10   ;   Ch = Ch_n10 
    546  
    547  
    548       !! II. Now starting iteration loop (IDM) 
    549       !! ------------------------------------- 
    550  
    551       DO jk = 1, nit 
    552   
    553 !CDIR NOVERRCHK 
    554          sqrt_Cd = SQRT( Cd ) 
    555  
    556          !! Turbulent scales : 
    557          !! ------------------ 
    558          U_star  = sqrt_Cd    * dU10             !  \\ L & Y eq. (7a) 
    559          T_star  = Ch/sqrt_Cd * dT               !  \\ L & Y eq. (7b) 
    560          q_star  = Ce/sqrt_Cd * dq               !  \\ L & Y eq. (7c) 
    561  
    562          !! Estimate the Monin-Obukov length : 
    563          !! ---------------------------------- 
    564          L  = (U_star*U_star) / ( vkarmn*grav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) ) 
    565  
    566          !! Stability parameters : 
    567          !! ---------------------- 
    568          zeta = zzu / L 
    569          zeta = SIGN( MIN( ABS( zeta ), 10.0 ), zeta ) 
    570  
    571          !! Psis, L & Y eq. (8c), (8d), (8e) : 
    572          !! ---------------------------------- 
    573 !CDIR NOVERRCHK 
    574          X2 = SQRT( ABS( 1. - 16.*zeta ) )   ;   X2 = MAX( X2 , 1.0 )    
    575 !CDIR NOVERRCHK 
    576          X  = SQRT( X2 ) 
    577  
    578          stabit    = 0.5 + SIGN( 0.5, zeta ) 
    579  
    580 !CDIR NOVERRCHK 
    581          psi_m = -5. * zeta * stabit   &                                                   ! Stable 
    582             &  + (1 - stabit)*(2*LOG((1. + X)/2) + LOG((1. + X2)/2) - 2*atan(X) + rpi/2)   ! Unstable 
    583  
    584 !CDIR NOVERRCHK 
    585          psi_h = -5. * zeta * stabit   &                                                   ! Stable 
    586             &  + (1 - stabit)*(2*LOG( (1. + X2)/2 ))                                       ! Unstable 
    587  
    588          !! Shifting the wind speed to 10m and neutral stability : 
    589          !! ------------------------------------------------------ 
    590 !CDIR NOVERRCHK 
    591          U_n10 = dU10 / (1. + SQRT( Cd_n10 ) / vkarmn * ( LOG(zzu/10.) - psi_m ) )    ! \\ L & Y eq. (9a) 
    592           
    593          !! Updating the neutral 10m transfer coefficients : 
    594          !! ------------------------------------------------ 
    595          Cd_n10  = 1e-3 * ( 2.7/U_n10 + 0.142 + U_n10/13.09 )              ! \\ L & Y eq. (6a) 
    596 !CDIR NOVERRCHK 
    597          sqrt_Cd_n10 = SQRT( Cd_n10 ) 
    598           
    599          Ce_n10  = 1e-3 * ( 34.6 * sqrt_Cd_n10 )                           ! \\ L & Y eq. (6b) 
    600          
    601          stab    = 0.5 + sign(0.5,zeta) 
    602          Ch_n10  = 1e-3 * sqrt_Cd_n10 * ( 18.*stab + 32.7*(1-stab) )       ! \\ L & Y eq. (6c), (6d) 
    603         
    604          !! Shifting the neutral  10m transfer coefficients to ( zzu , zeta ) : 
    605          !! -------------------------------------------------------------------- 
    606          !! Problem here, formulation used within L & Y differs from the one provided  
    607          !! in their fortran code (only for Ce and Ch) 
    608           
    609 !CDIR NOVERRCHK 
    610          Cd = Cd_n10/(1. + sqrt_Cd_n10/vkarmn*(LOG(zzu/10) - psi_m))**2     ! \\ L & Y eq. (10a) 
    611          
    612 !CDIR NOVERRCHK 
    613          xlogt = LOG(zzu/10) - psi_h 
    614 !CDIR NOVERRCHK 
    615          !?      Ch = Ch_n10*SQRT(Cd/Cd_n10)/(1. + Ch_n10/(vkarmn*sqrt_Cd_n10)*xlogt) 
    616          Ch = Ch_n10/( 1. + Ch_n10*xlogt/vkarmn/sqrt_Cd_n10 )**2             ! \\ L & Y eq. (10b) 
    617         
    618 !CDIR NOVERRCHK 
    619          !?      Ce = Ce_n10*SQRT(Cd/Cd_n10)/(1. + Ce_n10/(vkarmn*sqrt_Cd_n10)*xlogt) 
    620          Ce = Ce_n10/( 1. + Ce_n10*xlogt/vkarmn/sqrt_Cd_n10 )**2             ! \\ L & Y eq. (10c) 
    621        
    622       END DO 
    623  
    624       C_d(:,:) = Cd(:,:) 
    625       C_h(:,:) = Ch(:,:) 
    626       C_e(:,:) = Ce(:,:) 
    627  
    628    END SUBROUTINE TURB_CORE 
    629614   
    630615   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.