Changeset 758
- Timestamp:
- 2007-12-07T15:44:32+01:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_001_SBC/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r710 r758 174 174 ! ! input fields at the current time-step 175 175 176 !!gmIF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN177 178 CALL blk_oce_core( sst_m, ssu_m, ssv_m ) ! set the ocean surface fluxes179 180 !!gmENDIF176 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 181 181 ! ! using CORE bulk formulea 182 182 END SUBROUTINE sbc_blk_core … … 271 271 ! & Cd(:,:), Ch(:,:), Ce(:,:) ) 272 272 !!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 ) 276 277 277 278 ! ... umasked Momentum : utau, vtau at U- and V_points, resp. … … 465 466 END SUBROUTINE blk_ice_core 466 467 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 467 613 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 fluxes471 !! 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 D474 !!475 !! Momentum, Latent and sensible heat exchange coefficients476 !!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, Grenoble481 !! brodeau@hmg.inpg.fr482 !!---------------------------------------------------------------------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 indices494 INTEGER , PARAMETER :: nit = 3 ! number of iterations495 INTEGER , DIMENSION(jpi,jpj) :: stab ! stability test integer496 INTEGER , DIMENSION(jpi,jpj) :: stabit ! stability within iterative loop497 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 coefficient502 Ce_n10, & ! 10m neutral latent coefficient503 Ch_n10, & ! 10m neutral sensible coefficient504 Cd, & ! drag coefficient505 Ce, & ! latent coefficient506 Ch, & ! sensible coefficient507 sqrt_Cd_n10, & ! root square of Cd_n10508 sqrt_Cd, & ! root square of Cd509 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 zzu515 X2, X, &516 psi_m, &517 psi_h, &518 U_n10, & ! neutral wind velocity at 10m [m]519 xlogt520 !!---------------------------------------------------------------------521 522 !! I. Preliminary stuffs523 !! --------------------524 525 ! ... Air/sea differences526 dU10 = MAX( 0.5, dU ) ! we do not want to fall under 0.5 m/s527 dT = T_a - T_0 ! assuming that T_a is allready the potential temp. at zzu528 dq = q_a - q_sat529 530 ! ... Virtual potential temperature531 T_vpot = T_a * ( 1. + 0.608 * q_a )532 533 ! ... Computing Neutral Drag Coefficient534 !CDIR NOVERRCHK535 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 = 0542 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_n10546 547 548 !! II. Now starting iteration loop (IDM)549 !! -------------------------------------550 551 DO jk = 1, nit552 553 !CDIR NOVERRCHK554 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 / L569 zeta = SIGN( MIN( ABS( zeta ), 10.0 ), zeta )570 571 !! Psis, L & Y eq. (8c), (8d), (8e) :572 !! ----------------------------------573 !CDIR NOVERRCHK574 X2 = SQRT( ABS( 1. - 16.*zeta ) ) ; X2 = MAX( X2 , 1.0 )575 !CDIR NOVERRCHK576 X = SQRT( X2 )577 578 stabit = 0.5 + SIGN( 0.5, zeta )579 580 !CDIR NOVERRCHK581 psi_m = -5. * zeta * stabit & ! Stable582 & + (1 - stabit)*(2*LOG((1. + X)/2) + LOG((1. + X2)/2) - 2*atan(X) + rpi/2) ! Unstable583 584 !CDIR NOVERRCHK585 psi_h = -5. * zeta * stabit & ! Stable586 & + (1 - stabit)*(2*LOG( (1. + X2)/2 )) ! Unstable587 588 !! Shifting the wind speed to 10m and neutral stability :589 !! ------------------------------------------------------590 !CDIR NOVERRCHK591 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 NOVERRCHK597 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 provided607 !! in their fortran code (only for Ce and Ch)608 609 !CDIR NOVERRCHK610 Cd = Cd_n10/(1. + sqrt_Cd_n10/vkarmn*(LOG(zzu/10) - psi_m))**2 ! \\ L & Y eq. (10a)611 612 !CDIR NOVERRCHK613 xlogt = LOG(zzu/10) - psi_h614 !CDIR NOVERRCHK615 !? 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 NOVERRCHK619 !? 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 DO623 624 C_d(:,:) = Cd(:,:)625 C_h(:,:) = Ch(:,:)626 C_e(:,:) = Ce(:,:)627 628 END SUBROUTINE TURB_CORE629 614 630 615 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
