- Timestamp:
- 2008-04-02T15:58:41+02:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_001_SBC/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r874 r875 57 57 REAL(wp), PARAMETER :: Cice = 1.63e-3 ! transfer coefficient over ice 58 58 59 LOGICAL :: ln_2m = .FALSE. !: logical flag for height of air temp. and hum 60 REAL(wp) :: alpha_precip=1. !: multiplication factor for precipitation 61 59 62 !! * Substitutions 60 63 # include "domzgr_substitute.h90" … … 107 110 TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read 108 111 TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " " 109 NAMELIST/namsbc_core/ cn_dir, sn_wndi, sn_wndj, sn_humi, sn_qsr, &110 & sn_qlw , sn_tair, sn_prec, sn_snow112 NAMELIST/namsbc_core/ cn_dir, ln_2m, alpha_precip, sn_wndi, sn_wndj, sn_humi, sn_qsr, & 113 & sn_qlw , sn_tair, sn_prec, sn_snow 111 114 !!--------------------------------------------------------------------- 112 115 … … 158 161 WRITE(numout,*) '~~~~~~~~~~~~ ' 159 162 WRITE(numout,*) ' namsbc_core Namelist' 163 WRITE(numout,*) ' ln_2m = ', ln_2m 164 WRITE(numout,*) ' alpha_precip = ', alpha_precip 160 165 WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)' 161 166 DO jf = 1, jpfld … … 167 172 & ' starting record: ', sf(jf)%nstrec 168 173 END DO 174 IF( ln_2m ) THEN 175 WRITE(numout,*) ' Calling TURB_CORE_2Z for bulk transfert coefficients' 176 ELSE 177 WRITE(numout,*) ' Calling TURB_CORE_1Z for bulk transfert coefficients' 178 ENDIF 179 WRITE(numout,*) 180 ! 169 181 ENDIF 170 182 ! … … 216 228 REAL(wp), DIMENSION(jpi,jpj) :: Ce ! tansfert coefficient for evaporation (Q_lat) 217 229 REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin 230 REAL(wp), DIMENSION(jpi,jpj) :: zt_zu ! air temperature at wind speed height 231 REAL(wp), DIMENSION(jpi,jpj) :: zq_zu ! air spec. hum. at wind speed height 218 232 !!--------------------------------------------------------------------- 219 233 … … 266 280 267 281 ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point : 268 !!gm bug? a the compiling phase, add a copy in temporary arrays... ==> check perf! 269 ! CALL TURB_CORE( 10., zst (:,:), sf(jp_tair)%fnow(:,:), & 270 ! & zqsatw(:,:), sf(jp_humi)%fnow(:,:), zwind_speed_t(:,:), & 271 ! & Cd(:,:), Ch(:,:), Ce(:,:) ) 272 !!gm end 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 ) 277 282 IF( ln_2m ) THEN 283 !! If air temp. and spec. hum. are given at different height (2m) than wind (10m) : 284 CALL TURB_CORE_2Z(2.,10., zst , sf(jp_tair)%fnow, & 285 & zqsatw, sf(jp_humi)%fnow, zwind_speed_t, & 286 & Cd, Ch, Ce, zt_zu, zq_zu ) 287 ELSE 288 !! If air temp. and spec. hum. are given at same height than wind (10m) : 289 !gm bug? at the compiling phase, add a copy in temporary arrays... ==> check perf 290 ! CALL TURB_CORE_1Z( 10., zst (:,:), sf(jp_tair)%fnow(:,:), & 291 ! & zqsatw(:,:), sf(jp_humi)%fnow(:,:), zwind_speed_t(:,:), & 292 ! & Cd(:,:), Ch(:,:), Ce(:,:) ) 293 !gm bug 294 CALL TURB_CORE_1Z( 10., zst , sf(jp_tair)%fnow, & 295 & zqsatw, sf(jp_humi)%fnow, zwind_speed_t, & 296 & Cd, Ch, Ce ) 297 ENDIF 298 278 299 ! ... utau, vtau at U- and V_points, resp. 279 300 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines … … 291 312 ! Turbulent fluxes over ocean 292 313 ! ----------------------------- 293 !CDIR COLLAPSE 294 zevap(:,:) = MAX( 0.e0, rhoa *Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:) ) * zwind_speed_t(:,:) ) ! Evaporation 295 !CDIR COLLAPSE 296 zqsb (:,:) = rhoa*cpa*Ch(:,:)*( zst (:,:) - sf(jp_tair)%fnow(:,:) ) * zwind_speed_t(:,:) ! Sensible Heat 314 IF( ln_2m ) THEN 315 ! Values of temp. and hum. adjusted to 10m must be used instead of 2m values 316 zevap(:,:) = MAX( 0.e0, rhoa *Ce(:,:)*( zqsatw(:,:) - zq_zu(:,:) ) * zwind_speed_t(:,:) ) ! Evaporation 317 zqsb (:,:) = rhoa*cpa*Ch(:,:)*( zst (:,:) - zt_zu(:,:) ) * zwind_speed_t(:,:) ! Sensible Heat 318 ELSE 319 !CDIR COLLAPSE 320 zevap(:,:) = MAX( 0.e0, rhoa *Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:) ) * zwind_speed_t(:,:) ) ! Evaporation 321 !CDIR COLLAPSE 322 zqsb (:,:) = rhoa*cpa*Ch(:,:)*( zst (:,:) - sf(jp_tair)%fnow(:,:) ) * zwind_speed_t(:,:) ! Sensible Heat 323 ENDIF 297 324 !CDIR COLLAPSE 298 325 zqla (:,:) = Lv * zevap(:,:) ! Latent Heat … … 305 332 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) ! Downward Non Solar flux 306 333 !CDIR COLLAPSE 307 emp (:,:) = zevap(:,:) - sf(jp_prec)%fnow(:,:) * tmask(:,:,1)308 !CDIR COLLAPSE 309 emps(:,:) = zevap(:,:) - sf(jp_prec)%fnow(:,:) * tmask(:,:,1)334 emp (:,:) = zevap(:,:) - sf(jp_prec)%fnow(:,:) * alpha_precip * tmask(:,:,1) 335 !CDIR COLLAPSE 336 emps(:,:) = zevap(:,:) - sf(jp_prec)%fnow(:,:) * alpha_precip * tmask(:,:,1) 310 337 ! 311 338 END SUBROUTINE blk_oce_core … … 459 486 460 487 !CDIR COLLAPSE 461 p_tpr(:,:) = sf(jp_prec)%fnow(:,:) ! total precipitation [kg/m2/s]462 !CDIR COLLAPSE 463 p_spr(:,:) = sf(jp_snow)%fnow(:,:) ! solid precipitation [kg/m2/s]488 p_tpr(:,:) = sf(jp_prec)%fnow(:,:) * alpha_precip ! total precipitation [kg/m2/s] 489 !CDIR COLLAPSE 490 p_spr(:,:) = sf(jp_snow)%fnow(:,:) * alpha_precip ! solid precipitation [kg/m2/s] 464 491 ! 465 492 END SUBROUTINE blk_ice_core 466 493 494 467 495 SUBROUTINE TURB_CORE_1Z(zu, sst, T_a, q_sat, q_a, & 468 496 & dU, Cd, Ch, Ce ) … … 588 616 END SUBROUTINE TURB_CORE_1Z 589 617 618 619 SUBROUTINE TURB_CORE_2Z(zt, zu, sst, T_zt, q_sat, q_zt, dU, Cd, Ch, Ce, T_zu, q_zu) 620 !!---------------------------------------------------------------------- 621 !! *** ROUTINE turb_core *** 622 !! 623 !! ** Purpose : Computes turbulent transfert coefficients of surface 624 !! fluxes according to Large & Yeager (2004). 625 !! 626 !! ** 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 627 !! Momentum, Latent and sensible heat exchange coefficients 628 !! Caution: this procedure should only be used in cases when air 629 !! temperature (T_air) and air specific humidity (q_air) are at 2m 630 !! whereas wind (dU) is at 10m. 631 !! 632 !! References : 633 !! Large & Yeager, 2004 : ??? 634 !! History : 635 !! 9.0 ! 06-12 (L. Brodeau) Original code for 2Z 636 !!---------------------------------------------------------------------- 637 !! * Arguments 638 REAL(wp), INTENT(in) :: & 639 zt, & ! height for T_zt and q_zt [m] 640 zu ! height for dU [m] 641 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: & 642 sst, & ! sea surface temperature [Kelvin] 643 T_zt, & ! potential air temperature [Kelvin] 644 q_sat, & ! sea surface specific humidity [kg/kg] 645 q_zt, & ! specific air humidity [kg/kg] 646 dU ! relative wind module |U(zu)-U(0)| [m/s] 647 REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: & 648 Cd, & ! transfer coefficient for momentum (tau) 649 Ch, & ! transfer coefficient for sensible heat (Q_sens) 650 Ce, & ! transfert coefficient for evaporation (Q_lat) 651 T_zu, & ! air temp. shifted at zu [K] 652 q_zu ! spec. hum. shifted at zu [kg/kg] 653 654 !! * Local declarations 655 REAL(wp), DIMENSION(jpi,jpj) :: & 656 dU10, & ! dU [m/s] 657 dT, & ! air/sea temperature differeence [K] 658 dq, & ! air/sea humidity difference [K] 659 Cd_n10, & ! 10m neutral drag coefficient 660 Ce_n10, & ! 10m neutral latent coefficient 661 Ch_n10, & ! 10m neutral sensible coefficient 662 sqrt_Cd_n10, & ! root square of Cd_n10 663 sqrt_Cd, & ! root square of Cd 664 T_vpot_u, & ! virtual potential temperature [K] 665 T_star, & ! turbulent scale of tem. fluct. 666 q_star, & ! turbulent humidity of temp. fluct. 667 U_star, & ! turb. scale of velocity fluct. 668 L, & ! Monin-Obukov length [m] 669 zeta_u, & ! stability parameter at height zu 670 zeta_t, & ! stability parameter at height zt 671 U_n10, & ! neutral wind velocity at 10m [m] 672 xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m 673 674 INTEGER :: j_itt 675 INTEGER, PARAMETER :: nb_itt = 3 ! number of itterations 676 INTEGER, DIMENSION(jpi,jpj) :: & 677 & stab ! 1st stability test integer 678 REAL(wp), PARAMETER :: & 679 grav = 9.8, & ! gravity 680 kappa = 0.4 ! von Karman's constant 681 682 !! * Start 683 684 !! Initial air/sea differences 685 dU10 = max(0.5, dU) ! we don't want to fall under 0.5 m/s 686 dT = T_zt - sst 687 dq = q_zt - q_sat 688 689 !! Neutral Drag Coefficient : 690 stab = 0.5 + sign(0.5,dT) ! stab = 1 if dT > 0 -> STABLE 691 Cd_n10 = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 692 sqrt_Cd_n10 = sqrt(Cd_n10) 693 Ce_n10 = 1E-3*( 34.6 * sqrt_Cd_n10 ) 694 Ch_n10 = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1 - stab)) 695 696 !! Initializing transf. coeff. with their first guess neutral equivalents : 697 Cd = Cd_n10 ; Ce = Ce_n10 ; Ch = Ch_n10 ; sqrt_Cd = sqrt(Cd) 698 699 !! Initializing z_u values with z_t values : 700 T_zu = T_zt ; q_zu = q_zt 701 702 !! * Now starting iteration loop 703 DO j_itt=1, nb_itt 704 dT = T_zu - sst ; dq = q_zu - q_sat ! Updating air/sea differences 705 T_vpot_u = T_zu*(1. + 0.608*q_zu) ! Updating virtual potential temperature at zu 706 U_star = sqrt_Cd*dU10 ! Updating turbulent scales : (L & Y eq. (7)) 707 T_star = Ch/sqrt_Cd*dT ! 708 q_star = Ce/sqrt_Cd*dq ! 709 !! 710 L = (U_star*U_star) & ! Estimate the Monin-Obukov length at height zu 711 & / (kappa*grav/T_vpot_u*(T_star*(1.+0.608*q_zu) + 0.608*T_zu*q_star)) 712 !! Stability parameters : 713 zeta_u = zu/L ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) 714 zeta_t = zt/L ; zeta_t = sign( min(abs(zeta_t),10.0), zeta_t ) 715 zpsi_hu = psi_h(zeta_u) 716 zpsi_ht = psi_h(zeta_t) 717 zpsi_m = psi_m(zeta_u) 718 !! 719 !! Shifting the wind speed to 10m and neutral stability : (L & Y eq.(9a)) 720 ! U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u))) 721 U_n10 = dU10/(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) 722 !! 723 !! Shifting temperature and humidity at zu : (L & Y eq. (9b-9c)) 724 ! T_zu = T_zt - T_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t)) 725 T_zu = T_zt - T_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht) 726 ! q_zu = q_zt - q_star/kappa*(log(zt/zu) + psi_h(zeta_u) - psi_h(zeta_t)) 727 q_zu = q_zt - q_star/kappa*(log(zt/zu) + zpsi_hu - zpsi_ht) 728 !! 729 !! q_zu cannot have a negative value : forcing 0 730 stab = 0.5 + sign(0.5,q_zu) ; q_zu = stab*q_zu 731 !! 732 !! Updating the neutral 10m transfer coefficients : 733 Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 734 sqrt_Cd_n10 = sqrt(Cd_n10) 735 Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) 736 stab = 0.5 + sign(0.5,zeta_u) 737 Ch_n10 = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d) 738 !! 739 !! 740 !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) : 741 ! xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u)) 742 xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m) 743 Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd) 744 !! 745 ! xlogt = log(zu/10.) - psi_h(zeta_u) 746 xlogt = log(zu/10.) - zpsi_hu 747 !! 748 xct = 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10 749 Ch = Ch_n10*sqrt_Cd/sqrt_Cd_n10/xct 750 !! 751 xct = 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10 752 Ce = Ce_n10*sqrt_Cd/sqrt_Cd_n10/xct 753 !! 754 !! 755 END DO 756 !! 757 END SUBROUTINE TURB_CORE_2Z 758 759 590 760 FUNCTION psi_m(zta) !! Psis, L & Y eq. (8c), (8d), (8e) 591 761 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta
Note: See TracChangeset
for help on using the changeset viewer.
