Changeset 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
- Timestamp:
- 2012-01-28T17:44:18+01:00 (12 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r2777 r3294 14 14 !! 3.2 ! 2009-04 (B. Lemaire) Introduce iom_put 15 15 !! 3.3 ! 2010-10 (S. Masson) add diurnal cycle 16 !! 3.4 ! 2011-11 (C. Harris) Fill arrays required by CICE 16 17 !!---------------------------------------------------------------------- 17 18 … … 32 33 USE in_out_manager ! I/O manager 33 34 USE lib_mpp ! distribued memory computing library 35 USE wrk_nemo ! work arrays 36 USE timing ! Timing 34 37 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 35 38 USE prtctl ! Print control 36 #if defined key_lim3 39 USE sbcwave,ONLY : cdn_wave !wave module 40 #if defined key_lim3 || defined key_cice 37 41 USE sbc_ice ! Surface boundary condition: ice fields 38 42 #endif … … 43 47 PUBLIC sbc_blk_core ! routine called in sbcmod module 44 48 PUBLIC blk_ice_core ! routine called in sbc_ice_lim module 49 PUBLIC turb_core_2z ! routine calles in sbcblk_mfs module 45 50 46 51 INTEGER , PARAMETER :: jpfld = 9 ! maximum number of files to read … … 182 187 ! ! surface ocean fluxes computed with CLIO bulk formulea 183 188 IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce_core( sf, sst_m, ssu_m, ssv_m ) 189 190 #if defined key_cice 191 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 192 qlw_ice(:,:,1) = sf(jp_qlw)%fnow(:,:,1) 193 qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 194 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 195 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 196 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac 197 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac 198 wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1) 199 wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1) 200 ENDIF 201 #endif 184 202 ! 185 203 END SUBROUTINE sbc_blk_core … … 207 225 !! ** Nota : sf has to be a dummy argument for AGRIF on NEC 208 226 !!--------------------------------------------------------------------- 209 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released210 USE wrk_nemo, ONLY: zwnd_i => wrk_2d_1 , zwnd_j => wrk_2d_2 ! wind speed components at T-point211 USE wrk_nemo, ONLY: zqsatw => wrk_2d_3 ! specific humidity at pst212 USE wrk_nemo, ONLY: zqlw => wrk_2d_4 , zqsb => wrk_2d_5 ! long wave and sensible heat fluxes213 USE wrk_nemo, ONLY: zqla => wrk_2d_6 , zevap => wrk_2d_7 ! latent heat fluxes and evaporation214 USE wrk_nemo, ONLY: Cd => wrk_2d_8 ! transfer coefficient for momentum (tau)215 USE wrk_nemo, ONLY: Ch => wrk_2d_9 ! transfer coefficient for sensible heat (Q_sens)216 USE wrk_nemo, ONLY: Ce => wrk_2d_10 ! transfer coefficient for evaporation (Q_lat)217 USE wrk_nemo, ONLY: zst => wrk_2d_11 ! surface temperature in Kelvin218 USE wrk_nemo, ONLY: zt_zu => wrk_2d_12 ! air temperature at wind speed height219 USE wrk_nemo, ONLY: zq_zu => wrk_2d_13 ! air spec. hum. at wind speed height220 !221 227 TYPE(fld), INTENT(in), DIMENSION(:) :: sf ! input data 222 228 REAL(wp) , INTENT(in), DIMENSION(:,:) :: pst ! surface temperature [Celcius] … … 226 232 INTEGER :: ji, jj ! dummy loop indices 227 233 REAL(wp) :: zcoef_qsatw, zztmp ! local variable 234 REAL(wp), DIMENSION(:,:), POINTER :: zwnd_i, zwnd_j ! wind speed components at T-point 235 REAL(wp), DIMENSION(:,:), POINTER :: zqsatw ! specific humidity at pst 236 REAL(wp), DIMENSION(:,:), POINTER :: zqlw, zqsb ! long wave and sensible heat fluxes 237 REAL(wp), DIMENSION(:,:), POINTER :: zqla, zevap ! latent heat fluxes and evaporation 238 REAL(wp), DIMENSION(:,:), POINTER :: Cd ! transfer coefficient for momentum (tau) 239 REAL(wp), DIMENSION(:,:), POINTER :: Ch ! transfer coefficient for sensible heat (Q_sens) 240 REAL(wp), DIMENSION(:,:), POINTER :: Ce ! tansfert coefficient for evaporation (Q_lat) 241 REAL(wp), DIMENSION(:,:), POINTER :: zst ! surface temperature in Kelvin 242 REAL(wp), DIMENSION(:,:), POINTER :: zt_zu ! air temperature at wind speed height 243 REAL(wp), DIMENSION(:,:), POINTER :: zq_zu ! air spec. hum. at wind speed height 228 244 !!--------------------------------------------------------------------- 229 230 IF( wrk_in_use(2, 1,2,3,4,5,6,7,8,9,10,11,12,13) ) THEN 231 CALL ctl_stop('blk_oce_core: requested workspace arrays unavailable') ; RETURN 232 ENDIF 245 ! 246 IF( nn_timing == 1 ) CALL timing_start('blk_oce_core') 247 ! 248 CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap ) 249 CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu ) 233 250 ! 234 251 ! local scalars ( place there for vector optimisation purposes) … … 380 397 ENDIF 381 398 ! 382 IF( wrk_not_released(2, 1,2,3,4,5,6,7,8,9,10,11,12,13) ) & 383 CALL ctl_stop('blk_oce_core: failed to release workspace arrays') 399 CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j, zqsatw, zqlw, zqsb, zqla, zevap ) 400 CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu ) 401 ! 402 IF( nn_timing == 1 ) CALL timing_stop('blk_oce_core') 384 403 ! 385 404 END SUBROUTINE blk_oce_core … … 403 422 !! caution : the net upward water flux has with mm/day unit 404 423 !!--------------------------------------------------------------------- 405 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released406 USE wrk_nemo, ONLY: z_wnds_t => wrk_2d_1 ! wind speed ( = | U10m - U_ice | ) at T-point407 USE wrk_nemo, ONLY: wrk_3d_4 , wrk_3d_5 , wrk_3d_6 , wrk_3d_7408 !!409 424 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pst ! ice surface temperature (>0, =rt0 over land) [Kelvin] 410 425 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pui ! ice surface velocity (i- and i- components [m/s] … … 434 449 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 435 450 !! 451 REAL(wp), DIMENSION(:,:) , POINTER :: z_wnds_t ! wind speed ( = | U10m - U_ice | ) at T-point 436 452 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice 437 453 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice … … 439 455 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice 440 456 !!--------------------------------------------------------------------- 457 ! 458 IF( nn_timing == 1 ) CALL timing_start('blk_ice_core') 459 ! 460 CALL wrk_alloc( jpi,jpj, z_wnds_t ) 461 CALL wrk_alloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 441 462 442 463 ijpl = pdim ! number of ice categories 443 444 ! Set-up access to workspace arrays445 IF( wrk_in_use(2, 1) .OR. wrk_in_use(3, 4,5,6,7) ) THEN446 CALL ctl_stop('blk_ice_core: requested workspace arrays unavailable') ; RETURN447 ELSE IF(ijpl > jpk) THEN448 CALL ctl_stop('blk_ice_core: no. of ice categories > jpk so wrk_nemo 3D workspaces cannot be used.')449 RETURN450 END IF451 ! Set-up pointers to sub-arrays of workspaces452 z_qlw => wrk_3d_4(:,:,1:ijpl)453 z_qsb => wrk_3d_5(:,:,1:ijpl)454 z_dqlw => wrk_3d_6(:,:,1:ijpl)455 z_dqsb => wrk_3d_7(:,:,1:ijpl)456 464 457 465 ! local scalars ( place there for vector optimisation purposes) … … 606 614 ENDIF 607 615 608 IF( wrk_not_released(2, 1) .OR. & 609 wrk_not_released(3, 4,5,6,7) ) CALL ctl_stop('blk_ice_core: failed to release workspace arrays') 616 CALL wrk_dealloc( jpi,jpj, z_wnds_t ) 617 CALL wrk_dealloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 618 ! 619 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core') 610 620 ! 611 621 END SUBROUTINE blk_ice_core … … 628 638 !! References : Large & Yeager, 2004 : ??? 629 639 !!---------------------------------------------------------------------- 630 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released631 USE wrk_nemo, ONLY: dU10 => wrk_2d_14 ! dU [m/s]632 USE wrk_nemo, ONLY: dT => wrk_2d_15 ! air/sea temperature difference [K]633 USE wrk_nemo, ONLY: dq => wrk_2d_16 ! air/sea humidity difference [K]634 USE wrk_nemo, ONLY: Cd_n10 => wrk_2d_17 ! 10m neutral drag coefficient635 USE wrk_nemo, ONLY: Ce_n10 => wrk_2d_18 ! 10m neutral latent coefficient636 USE wrk_nemo, ONLY: Ch_n10 => wrk_2d_19 ! 10m neutral sensible coefficient637 USE wrk_nemo, ONLY: sqrt_Cd_n10 => wrk_2d_20 ! root square of Cd_n10638 USE wrk_nemo, ONLY: sqrt_Cd => wrk_2d_21 ! root square of Cd639 USE wrk_nemo, ONLY: T_vpot => wrk_2d_22 ! virtual potential temperature [K]640 USE wrk_nemo, ONLY: T_star => wrk_2d_23 ! turbulent scale of tem. fluct.641 USE wrk_nemo, ONLY: q_star => wrk_2d_24 ! turbulent humidity of temp. fluct.642 USE wrk_nemo, ONLY: U_star => wrk_2d_25 ! turb. scale of velocity fluct.643 USE wrk_nemo, ONLY: L => wrk_2d_26 ! Monin-Obukov length [m]644 USE wrk_nemo, ONLY: zeta => wrk_2d_27 ! stability parameter at height zu645 USE wrk_nemo, ONLY: U_n10 => wrk_2d_28 ! neutral wind velocity at 10m [m]646 USE wrk_nemo, ONLY: xlogt => wrk_2d_29, xct => wrk_2d_30, &647 zpsi_h => wrk_2d_31, zpsi_m => wrk_2d_32648 USE wrk_nemo, ONLY: stab => iwrk_2d_1 ! 1st guess stability test integer649 !650 640 REAL(wp) , INTENT(in ) :: zu ! altitude of wind measurement [m] 651 641 REAL(wp), DIMENSION(:,:), INTENT(in ) :: sst ! sea surface temperature [Kelvin] … … 662 652 REAL(wp), PARAMETER :: grav = 9.8 ! gravity 663 653 REAL(wp), PARAMETER :: kappa = 0.4 ! von Karman s constant 654 655 REAL(wp), DIMENSION(:,:), POINTER :: dU10 ! dU [m/s] 656 REAL(wp), DIMENSION(:,:), POINTER :: dT ! air/sea temperature differeence [K] 657 REAL(wp), DIMENSION(:,:), POINTER :: dq ! air/sea humidity difference [K] 658 REAL(wp), DIMENSION(:,:), POINTER :: Cd_n10 ! 10m neutral drag coefficient 659 REAL(wp), DIMENSION(:,:), POINTER :: Ce_n10 ! 10m neutral latent coefficient 660 REAL(wp), DIMENSION(:,:), POINTER :: Ch_n10 ! 10m neutral sensible coefficient 661 REAL(wp), DIMENSION(:,:), POINTER :: sqrt_Cd_n10 ! root square of Cd_n10 662 REAL(wp), DIMENSION(:,:), POINTER :: sqrt_Cd ! root square of Cd 663 REAL(wp), DIMENSION(:,:), POINTER :: T_vpot ! virtual potential temperature [K] 664 REAL(wp), DIMENSION(:,:), POINTER :: T_star ! turbulent scale of tem. fluct. 665 REAL(wp), DIMENSION(:,:), POINTER :: q_star ! turbulent humidity of temp. fluct. 666 REAL(wp), DIMENSION(:,:), POINTER :: U_star ! turb. scale of velocity fluct. 667 REAL(wp), DIMENSION(:,:), POINTER :: L ! Monin-Obukov length [m] 668 REAL(wp), DIMENSION(:,:), POINTER :: zeta ! stability parameter at height zu 669 REAL(wp), DIMENSION(:,:), POINTER :: U_n10 ! neutral wind velocity at 10m [m] 670 REAL(wp), DIMENSION(:,:), POINTER :: xlogt, xct, zpsi_h, zpsi_m 671 672 INTEGER , DIMENSION(:,:), POINTER :: stab ! 1st guess stability test integer 664 673 !!---------------------------------------------------------------------- 665 666 IF( wrk_in_use(2, 14,15,16,17,18,19, & 667 20,21,22,23,24,25,26,27,28,29, & 668 30,31,32) .OR. & 669 iwrk_in_use(2, 1) ) THEN 670 CALL ctl_stop('TURB_CORE_1Z: requested workspace arrays unavailable') ; RETURN 671 ENDIF 674 ! 675 IF( nn_timing == 1 ) CALL timing_start('TURB_CORE_1Z') 676 ! 677 CALL wrk_alloc( jpi,jpj, stab ) ! integer 678 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 679 CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 672 680 673 681 !! * Start … … 682 690 !! Neutral Drag Coefficient 683 691 stab = 0.5 + sign(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0 684 Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 692 IF ( ln_cdgw ) THEN 693 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 694 Cd_n10(:,:) = cdn_wave 695 ELSE 696 Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! L & Y eq. (6a) 697 ENDIF 685 698 sqrt_Cd_n10 = sqrt(Cd_n10) 686 699 Ce_n10 = 1E-3 * ( 34.6 * sqrt_Cd_n10 ) ! L & Y eq. (6b) … … 705 718 zpsi_m = psi_m(zeta) 706 719 707 !! Shifting the wind speed to 10m and neutral stability : 708 U_n10 = dU10*1./(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) ! L & Y eq. (9a) 709 710 !! Updating the neutral 10m transfer coefficients : 711 Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 712 sqrt_Cd_n10 = sqrt(Cd_n10) 713 Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) 714 stab = 0.5 + sign(0.5,zeta) 715 Ch_n10 = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c), (6d) 716 717 !! Shifting the neutral 10m transfer coefficients to ( zu , zeta ) : 718 !! 719 xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10) - zpsi_m) 720 Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd) 720 IF ( ln_cdgw ) THEN 721 sqrt_Cd=kappa/((kappa/sqrt_Cd_n10) - zpsi_m) ; Cd=sqrt_Cd*sqrt_Cd; 722 ELSE 723 !! Shifting the wind speed to 10m and neutral stability : 724 U_n10 = dU10*1./(1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m)) ! L & Y eq. (9a) 725 726 !! Updating the neutral 10m transfer coefficients : 727 Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 728 sqrt_Cd_n10 = sqrt(Cd_n10) 729 Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) 730 stab = 0.5 + sign(0.5,zeta) 731 Ch_n10 = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c), (6d) 732 733 !! Shifting the neutral 10m transfer coefficients to ( zu , zeta ) : 734 !! 735 xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10) - zpsi_m) 736 Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd) 737 ENDIF 721 738 !! 722 739 xlogt = log(zu/10.) - zpsi_h … … 730 747 END DO 731 748 !! 732 IF( wrk_not_released(2, 14,15,16,17,18,19, &733 & 20,21,22,23,24,25,26,27,28,29, &734 & 30,31,32 ) .OR. &735 iwrk_not_released(2, 1) ) &736 CALL ctl_stop('TURB_CORE_1Z: failed to release workspace arrays')749 CALL wrk_dealloc( jpi,jpj, stab ) ! integer 750 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 751 CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 752 ! 753 IF( nn_timing == 1 ) CALL timing_stop('TURB_CORE_1Z') 737 754 ! 738 755 END SUBROUTINE TURB_CORE_1Z … … 754 771 !! References : Large & Yeager, 2004 : ??? 755 772 !!---------------------------------------------------------------------- 756 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released 757 USE wrk_nemo, ONLY: dU10 => wrk_2d_14 ! dU [m/s] 758 USE wrk_nemo, ONLY: dT => wrk_2d_15 ! air/sea temperature difference [K] 759 USE wrk_nemo, ONLY: dq => wrk_2d_16 ! air/sea humidity difference [K] 760 USE wrk_nemo, ONLY: Cd_n10 => wrk_2d_17 ! 10m neutral drag coefficient 761 USE wrk_nemo, ONLY: Ce_n10 => wrk_2d_18 ! 10m neutral latent coefficient 762 USE wrk_nemo, ONLY: Ch_n10 => wrk_2d_19 ! 10m neutral sensible coefficient 763 USE wrk_nemo, ONLY: sqrt_Cd_n10 => wrk_2d_20 ! root square of Cd_n10 764 USE wrk_nemo, ONLY: sqrt_Cd => wrk_2d_21 ! root square of Cd 765 USE wrk_nemo, ONLY: T_vpot => wrk_2d_22 ! virtual potential temperature [K] 766 USE wrk_nemo, ONLY: T_star => wrk_2d_23 ! turbulent scale of tem. fluct. 767 USE wrk_nemo, ONLY: q_star => wrk_2d_24 ! turbulent humidity of temp. fluct. 768 USE wrk_nemo, ONLY: U_star => wrk_2d_25 ! turb. scale of velocity fluct. 769 USE wrk_nemo, ONLY: L => wrk_2d_26 ! Monin-Obukov length [m] 770 USE wrk_nemo, ONLY: zeta_u => wrk_2d_27 ! stability parameter at height zu 771 USE wrk_nemo, ONLY: zeta_t => wrk_2d_28 ! stability parameter at height zt 772 USE wrk_nemo, ONLY: U_n10 => wrk_2d_29 ! neutral wind velocity at 10m [m] 773 USE wrk_nemo, ONLY: xlogt => wrk_2d_30, xct => wrk_2d_31, zpsi_hu => wrk_2d_32, zpsi_ht => wrk_2d_33, zpsi_m => wrk_2d_34 774 USE wrk_nemo, ONLY: stab => iwrk_2d_1 ! 1st guess stability test integer 775 !! 776 REAL(wp), INTENT(in) :: & 777 zt, & ! height for T_zt and q_zt [m] 778 zu ! height for dU [m] 779 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: & 780 sst, & ! sea surface temperature [Kelvin] 781 T_zt, & ! potential air temperature [Kelvin] 782 q_sat, & ! sea surface specific humidity [kg/kg] 783 q_zt, & ! specific air humidity [kg/kg] 784 dU ! relative wind module |U(zu)-U(0)| [m/s] 785 REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: & 786 Cd, & ! transfer coefficient for momentum (tau) 787 Ch, & ! transfer coefficient for sensible heat (Q_sens) 788 Ce, & ! transfert coefficient for evaporation (Q_lat) 789 T_zu, & ! air temp. shifted at zu [K] 790 q_zu ! spec. hum. shifted at zu [kg/kg] 773 REAL(wp), INTENT(in ) :: zt ! height for T_zt and q_zt [m] 774 REAL(wp), INTENT(in ) :: zu ! height for dU [m] 775 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst ! sea surface temperature [Kelvin] 776 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: T_zt ! potential air temperature [Kelvin] 777 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_sat ! sea surface specific humidity [kg/kg] 778 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity [kg/kg] 779 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: dU ! relative wind module |U(zu)-U(0)| [m/s] 780 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) 781 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens) 782 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Ce ! transfert coefficient for evaporation (Q_lat) 783 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: T_zu ! air temp. shifted at zu [K] 784 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: q_zu ! spec. hum. shifted at zu [kg/kg] 791 785 792 786 INTEGER :: j_itt 793 INTEGER, PARAMETER :: nb_itt = 3 ! number of itterations 794 REAL(wp), PARAMETER :: & 795 grav = 9.8, & ! gravity 796 kappa = 0.4 ! von Karman's constant 787 INTEGER , PARAMETER :: nb_itt = 3 ! number of itterations 788 REAL(wp), PARAMETER :: grav = 9.8 ! gravity 789 REAL(wp), PARAMETER :: kappa = 0.4 ! von Karman's constant 790 791 REAL(wp), DIMENSION(:,:), POINTER :: dU10 ! dU [m/s] 792 REAL(wp), DIMENSION(:,:), POINTER :: dT ! air/sea temperature differeence [K] 793 REAL(wp), DIMENSION(:,:), POINTER :: dq ! air/sea humidity difference [K] 794 REAL(wp), DIMENSION(:,:), POINTER :: Cd_n10 ! 10m neutral drag coefficient 795 REAL(wp), DIMENSION(:,:), POINTER :: Ce_n10 ! 10m neutral latent coefficient 796 REAL(wp), DIMENSION(:,:), POINTER :: Ch_n10 ! 10m neutral sensible coefficient 797 REAL(wp), DIMENSION(:,:), POINTER :: sqrt_Cd_n10 ! root square of Cd_n10 798 REAL(wp), DIMENSION(:,:), POINTER :: sqrt_Cd ! root square of Cd 799 REAL(wp), DIMENSION(:,:), POINTER :: T_vpot ! virtual potential temperature [K] 800 REAL(wp), DIMENSION(:,:), POINTER :: T_star ! turbulent scale of tem. fluct. 801 REAL(wp), DIMENSION(:,:), POINTER :: q_star ! turbulent humidity of temp. fluct. 802 REAL(wp), DIMENSION(:,:), POINTER :: U_star ! turb. scale of velocity fluct. 803 REAL(wp), DIMENSION(:,:), POINTER :: L ! Monin-Obukov length [m] 804 REAL(wp), DIMENSION(:,:), POINTER :: zeta_u ! stability parameter at height zu 805 REAL(wp), DIMENSION(:,:), POINTER :: zeta_t ! stability parameter at height zt 806 REAL(wp), DIMENSION(:,:), POINTER :: U_n10 ! neutral wind velocity at 10m [m] 807 REAL(wp), DIMENSION(:,:), POINTER :: xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m 808 809 INTEGER , DIMENSION(:,:), POINTER :: stab ! 1st stability test integer 797 810 !!---------------------------------------------------------------------- 798 !! * Start 799 800 IF( wrk_in_use(2, 14,15,16,17,18,19, & 801 20,21,22,23,24,25,26,27,28,29, & 802 30,31,32,33,34) .OR. & 803 iwrk_in_use(2, 1) ) THEN 804 CALL ctl_stop('TURB_CORE_2Z: requested workspace arrays unavailable') ; RETURN 805 ENDIF 811 ! 812 IF( nn_timing == 1 ) CALL timing_start('TURB_CORE_2Z') 813 ! 814 CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 815 CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 816 CALL wrk_alloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m ) 817 CALL wrk_alloc( jpi,jpj, stab ) ! interger 806 818 807 819 !! Initial air/sea differences … … 812 824 !! Neutral Drag Coefficient : 813 825 stab = 0.5 + sign(0.5,dT) ! stab = 1 if dT > 0 -> STABLE 814 Cd_n10 = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 826 IF( ln_cdgw ) THEN 827 cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 828 Cd_n10(:,:) = cdn_wave 829 ELSE 830 Cd_n10 = 1E-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 831 ENDIF 815 832 sqrt_Cd_n10 = sqrt(Cd_n10) 816 833 Ce_n10 = 1E-3*( 34.6 * sqrt_Cd_n10 ) … … 853 870 stab = 0.5 + sign(0.5,q_zu) ; q_zu = stab*q_zu 854 871 !! 855 !! Updating the neutral 10m transfer coefficients : 856 Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 857 sqrt_Cd_n10 = sqrt(Cd_n10) 858 Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) 859 stab = 0.5 + sign(0.5,zeta_u) 860 Ch_n10 = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d) 861 !! 862 !! 863 !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) : 864 ! xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - psi_m(zeta_u)) 865 xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m) 866 Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd) 867 !! 868 ! xlogt = log(zu/10.) - psi_h(zeta_u) 872 IF( ln_cdgw ) THEN 873 sqrt_Cd=kappa/((kappa/sqrt_Cd_n10) - zpsi_m) ; Cd=sqrt_Cd*sqrt_Cd; 874 ELSE 875 !! Updating the neutral 10m transfer coefficients : 876 Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! L & Y eq. (6a) 877 sqrt_Cd_n10 = sqrt(Cd_n10) 878 Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! L & Y eq. (6b) 879 stab = 0.5 + sign(0.5,zeta_u) 880 Ch_n10 = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! L & Y eq. (6c-6d) 881 !! 882 !! 883 !! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) : 884 xct = 1. + sqrt_Cd_n10/kappa*(log(zu/10.) - zpsi_m) 885 Cd = Cd_n10/(xct*xct) ; sqrt_Cd = sqrt(Cd) 886 ENDIF 887 !! 869 888 xlogt = log(zu/10.) - zpsi_hu 870 889 !! … … 878 897 END DO 879 898 !! 880 IF( wrk_not_released(2, 14,15,16,17,18,19, & 881 & 20,21,22,23,24,25,26,27,28,29, & 882 & 30,31,32,33,34 ) .OR. & 883 iwrk_not_released(2, 1) ) & 884 CALL ctl_stop('TURB_CORE_2Z: failed to release workspace arrays') 899 CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 900 CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 901 CALL wrk_dealloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m ) 902 CALL wrk_dealloc( jpi,jpj, stab ) ! interger 903 ! 904 IF( nn_timing == 1 ) CALL timing_stop('TURB_CORE_2Z') 885 905 ! 886 906 END SUBROUTINE TURB_CORE_2Z … … 889 909 FUNCTION psi_m(zta) !! Psis, L & Y eq. (8c), (8d), (8e) 890 910 !------------------------------------------------------------------------------- 891 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released892 USE wrk_nemo, ONLY: X2 => wrk_2d_35893 USE wrk_nemo, ONLY: X => wrk_2d_36894 USE wrk_nemo, ONLY: stabit => wrk_2d_37895 !!896 911 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 897 912 898 913 REAL(wp), PARAMETER :: pi = 3.141592653589793_wp 899 914 REAL(wp), DIMENSION(jpi,jpj) :: psi_m 915 REAL(wp), DIMENSION(:,:), POINTER :: X2, X, stabit 900 916 !------------------------------------------------------------------------------- 901 917 902 IF( wrk_in_use(2, 35,36,37) ) THEN 903 CALL ctl_stop('psi_m: requested workspace arrays unavailable') ; RETURN 904 ENDIF 918 CALL wrk_alloc( jpi,jpj, X2, X, stabit ) 905 919 906 920 X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.0) ; X = sqrt(X2) … … 909 923 & + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable 910 924 911 IF( wrk_not_released(2, 35,36,37) ) CALL ctl_stop('psi_m: failed to release workspace arrays')925 CALL wrk_dealloc( jpi,jpj, X2, X, stabit ) 912 926 ! 913 927 END FUNCTION psi_m … … 916 930 FUNCTION psi_h( zta ) !! Psis, L & Y eq. (8c), (8d), (8e) 917 931 !------------------------------------------------------------------------------- 918 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released919 USE wrk_nemo, ONLY: X2 => wrk_2d_35920 USE wrk_nemo, ONLY: X => wrk_2d_36921 USE wrk_nemo, ONLY: stabit => wrk_2d_37922 !923 932 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 924 933 ! 925 934 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 935 REAL(wp), DIMENSION(:,:), POINTER :: X2, X, stabit 926 936 !------------------------------------------------------------------------------- 927 937 928 IF( wrk_in_use(2, 35,36,37) ) THEN 929 CALL ctl_stop('psi_h: requested workspace arrays unavailable') ; RETURN 930 ENDIF 938 CALL wrk_alloc( jpi,jpj, X2, X, stabit ) 931 939 932 940 X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.) ; X = sqrt(X2) … … 935 943 & + (1. - stabit)*(2.*log( (1. + X2)/2. )) ! Unstable 936 944 937 IF( wrk_not_released(2, 35,36,37) ) CALL ctl_stop('psi_h: failed to release workspace arrays')945 CALL wrk_dealloc( jpi,jpj, X2, X, stabit ) 938 946 ! 939 947 END FUNCTION psi_h
Note: See TracChangeset
for help on using the changeset viewer.
