Changeset 2715 for trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
- Timestamp:
- 2011-03-30T17:58:35+02:00 (13 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r2528 r2715 11 11 !! - Implement reading of 6-hourly fields 12 12 !! 3.0 ! 2006-06 (G. Madec) sbc rewritting 13 !! - ! 2006-12 (L. Brodeau) Original code for TURB_CORE_2Z 13 14 !! 3.2 ! 2009-04 (B. Lemaire) Introduce iom_put 14 15 !! 3.3 ! 2010-10 (S. Masson) add diurnal cycle … … 40 41 PRIVATE 41 42 42 PUBLIC sbc_blk_core ! routine called in sbcmod module43 PUBLIC blk_ice_core ! routine called in sbc_ice_lim module44 43 PUBLIC sbc_blk_core ! routine called in sbcmod module 44 PUBLIC blk_ice_core ! routine called in sbc_ice_lim module 45 45 46 INTEGER , PARAMETER :: jpfld = 9 ! maximum number of files to read 46 47 INTEGER , PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point … … 53 54 INTEGER , PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s) 54 55 INTEGER , PARAMETER :: jp_tdif = 9 ! index of tau diff associated to HF tau (N/m2) at T-point 56 55 57 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read) 56 58 … … 110 112 !! - emp, emps evaporation minus precipitation 111 113 !!---------------------------------------------------------------------- 112 INTEGER, INTENT( in) :: kt ! ocean time step114 INTEGER, INTENT(in) :: kt ! ocean time step 113 115 !! 114 116 INTEGER :: ierror ! return error code … … 166 168 ! 167 169 ALLOCATE( sf(jfld), STAT=ierror ) ! set sf structure 168 IF( ierror > 0 ) THEN 169 CALL ctl_stop( 'sbc_blk_core: unable to allocate sf structure' ) ; RETURN 170 ENDIF 170 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_core: unable to allocate sf structure' ) 171 171 DO ifpr= 1, jfld 172 172 ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 173 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )173 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 174 174 END DO 175 175 ! ! fill sf with slf_i and control print … … 210 210 !! ** Nota : sf has to be a dummy argument for AGRIF on NEC 211 211 !!--------------------------------------------------------------------- 212 TYPE(fld), INTENT(in), DIMENSION(:) :: sf ! input data 213 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: pst ! surface temperature [Celcius] 214 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: pu ! surface current at U-point (i-component) [m/s] 215 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: pv ! surface current at V-point (j-component) [m/s] 216 217 INTEGER :: ji, jj ! dummy loop indices 218 REAL(wp) :: zcoef_qsatw 219 REAL(wp) :: zztmp ! temporary variable 220 REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point 221 REAL(wp), DIMENSION(jpi,jpj) :: zqsatw ! specific humidity at pst 222 REAL(wp), DIMENSION(jpi,jpj) :: zqlw, zqsb ! long wave and sensible heat fluxes 223 REAL(wp), DIMENSION(jpi,jpj) :: zqla, zevap ! latent heat fluxes and evaporation 224 REAL(wp), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) 225 REAL(wp), DIMENSION(jpi,jpj) :: Ch ! transfer coefficient for sensible heat (Q_sens) 226 REAL(wp), DIMENSION(jpi,jpj) :: Ce ! tansfert coefficient for evaporation (Q_lat) 227 REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin 228 REAL(wp), DIMENSION(jpi,jpj) :: zt_zu ! air temperature at wind speed height 229 REAL(wp), DIMENSION(jpi,jpj) :: zq_zu ! air spec. hum. at wind speed height 212 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 213 USE wrk_nemo, ONLY: zwnd_i => wrk_2d_1 , zwnd_j => wrk_2d_2 ! wind speed components at T-point 214 USE wrk_nemo, ONLY: zqsatw => wrk_2d_3 ! specific humidity at pst 215 USE wrk_nemo, ONLY: zqlw => wrk_2d_4 , zqsb => wrk_2d_5 ! long wave and sensible heat fluxes 216 USE wrk_nemo, ONLY: zqla => wrk_2d_6 , zevap => wrk_2d_7 ! latent heat fluxes and evaporation 217 USE wrk_nemo, ONLY: Cd => wrk_2d_8 ! transfer coefficient for momentum (tau) 218 USE wrk_nemo, ONLY: Ch => wrk_2d_9 ! transfer coefficient for sensible heat (Q_sens) 219 USE wrk_nemo, ONLY: Ce => wrk_2d_10 ! transfer coefficient for evaporation (Q_lat) 220 USE wrk_nemo, ONLY: zst => wrk_2d_11 ! surface temperature in Kelvin 221 USE wrk_nemo, ONLY: zt_zu => wrk_2d_12 ! air temperature at wind speed height 222 USE wrk_nemo, ONLY: zq_zu => wrk_2d_13 ! air spec. hum. at wind speed height 223 ! 224 TYPE(fld), INTENT(in), DIMENSION(:) :: sf ! input data 225 REAL(wp) , INTENT(in), DIMENSION(:,:) :: pst ! surface temperature [Celcius] 226 REAL(wp) , INTENT(in), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 227 REAL(wp) , INTENT(in), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 228 ! 229 INTEGER :: ji, jj ! dummy loop indices 230 REAL(wp) :: zcoef_qsatw, zztmp ! local variable 230 231 !!--------------------------------------------------------------------- 231 232 233 IF( wrk_in_use(2, 1,2,3,4,5,6,7,8,9,10,11,12,13) ) THEN 234 CALL ctl_stop('blk_oce_core: requested workspace arrays unavailable') ; RETURN 235 ENDIF 236 ! 232 237 ! local scalars ( place there for vector optimisation purposes) 233 238 zcoef_qsatw = 0.98 * 640380. / rhoa … … 293 298 ! & Cd (:,:), Ch (:,:), Ce (:,:) ) 294 299 !gm bug 295 CALL TURB_CORE_1Z( 10., zst , sf(jp_tair)%fnow, & 296 & zqsatw, sf(jp_humi)%fnow, wndm, & 300 ! ARPDBG - this won't compile with gfortran. Fix but check performance 301 ! as per comment above. 302 CALL TURB_CORE_1Z( 10., zst , sf(jp_tair)%fnow(:,:,1), & 303 & zqsatw, sf(jp_humi)%fnow(:,:,1), wndm, & 297 304 & Cd , Ch , Ce ) 298 305 ENDIF … … 376 383 ENDIF 377 384 ! 385 IF( wrk_not_released(2, 1,2,3,4,5,6,7,8,9,10,11,12,13) ) & 386 CALL ctl_stop('blk_oce_core: failed to release workspace arrays') 387 ! 378 388 END SUBROUTINE blk_oce_core 379 389 … … 396 406 !! caution : the net upward water flux has with mm/day unit 397 407 !!--------------------------------------------------------------------- 398 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pst ! ice surface temperature (>0, =rt0 over land) [Kelvin] 399 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pui ! ice surface velocity (i- and i- components [m/s] 400 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvi ! at I-point (B-grid) or U & V-point (C-grid) 401 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: palb ! ice albedo (clear sky) (alb_ice_cs) [%] 402 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_taui ! i- & j-components of surface ice stress [N/m2] 403 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_tauj ! at I-point (B-grid) or U & V-point (C-grid) 404 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: p_qns ! non solar heat flux over ice (T-point) [W/m2] 405 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: p_qsr ! solar heat flux over ice (T-point) [W/m2] 406 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: p_qla ! latent heat flux over ice (T-point) [W/m2] 407 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: p_dqns ! non solar heat sensistivity (T-point) [W/m2] 408 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: p_dqla ! latent heat sensistivity (T-point) [W/m2] 409 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_tpr ! total precipitation (T-point) [Kg/m2/s] 410 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_spr ! solid precipitation (T-point) [Kg/m2/s] 411 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_fr1 ! 1sr fraction of qsr penetration in ice (T-point) [%] 412 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: p_fr2 ! 2nd fraction of qsr penetration in ice (T-point) [%] 413 CHARACTER(len=1) , INTENT(in ) :: cd_grid ! ice grid ( C or B-grid) 414 INTEGER , INTENT(in ) :: pdim ! number of ice categories 408 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 409 USE wrk_nemo, ONLY: z_wnds_t => wrk_2d_1 ! wind speed ( = | U10m - U_ice | ) at T-point 410 USE wrk_nemo, ONLY: wrk_3d_4 , wrk_3d_5 , wrk_3d_6 , wrk_3d_7 411 !! 412 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pst ! ice surface temperature (>0, =rt0 over land) [Kelvin] 413 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pui ! ice surface velocity (i- and i- components [m/s] 414 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvi ! at I-point (B-grid) or U & V-point (C-grid) 415 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb ! ice albedo (clear sky) (alb_ice_cs) [%] 416 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_taui ! i- & j-components of surface ice stress [N/m2] 417 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_tauj ! at I-point (B-grid) or U & V-point (C-grid) 418 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qns ! non solar heat flux over ice (T-point) [W/m2] 419 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qsr ! solar heat flux over ice (T-point) [W/m2] 420 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qla ! latent heat flux over ice (T-point) [W/m2] 421 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_dqns ! non solar heat sensistivity (T-point) [W/m2] 422 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_dqla ! latent heat sensistivity (T-point) [W/m2] 423 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_tpr ! total precipitation (T-point) [Kg/m2/s] 424 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_spr ! solid precipitation (T-point) [Kg/m2/s] 425 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_fr1 ! 1sr fraction of qsr penetration in ice (T-point) [%] 426 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_fr2 ! 2nd fraction of qsr penetration in ice (T-point) [%] 427 CHARACTER(len=1) , INTENT(in ) :: cd_grid ! ice grid ( C or B-grid) 428 INTEGER , INTENT(in ) :: pdim ! number of ice categories 415 429 !! 416 430 INTEGER :: ji, jj, jl ! dummy loop indices … … 422 436 REAL(wp) :: zwnorm_f, zwndi_f , zwndj_f ! relative wind module and components at F-point 423 437 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 424 REAL(wp), DIMENSION(jpi,jpj) :: z_wnds_t ! wind speed ( = | U10m - U_ice | ) at T-point425 REAL(wp), DIMENSION( jpi,jpj,pdim) :: z_qlw! long wave heat flux over ice426 REAL(wp), DIMENSION( jpi,jpj,pdim) :: z_qsb! sensible heat flux over ice427 REAL(wp), DIMENSION( jpi,jpj,pdim) :: z_dqlw! long wave heat sensitivity over ice428 REAL(wp), DIMENSION( jpi,jpj,pdim) :: z_dqsb! sensible heat sensitivity over ice438 !! 439 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice 440 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice 441 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqlw ! long wave heat sensitivity over ice 442 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice 429 443 !!--------------------------------------------------------------------- 430 444 431 445 ijpl = pdim ! number of ice categories 446 447 ! Set-up access to workspace arrays 448 IF( wrk_in_use(2, 1) .OR. wrk_in_use(3, 4,5,6,7) ) THEN 449 CALL ctl_stop('blk_ice_core: requested workspace arrays unavailable') ; RETURN 450 ELSE IF(ijpl > jpk) THEN 451 CALL ctl_stop('blk_ice_core: no. of ice categories > jpk so wrk_nemo 3D workspaces cannot be used.') 452 RETURN 453 END IF 454 ! Set-up pointers to sub-arrays of workspaces 455 z_qlw => wrk_3d_4(:,:,1:ijpl) 456 z_qsb => wrk_3d_5(:,:,1:ijpl) 457 z_dqlw => wrk_3d_6(:,:,1:ijpl) 458 z_dqsb => wrk_3d_7(:,:,1:ijpl) 432 459 433 460 ! local scalars ( place there for vector optimisation purposes) … … 579 606 ENDIF 580 607 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') 610 ! 581 611 END SUBROUTINE blk_ice_core 582 612 583 613 584 614 SUBROUTINE TURB_CORE_1Z(zu, sst, T_a, q_sat, q_a, & 585 & dU , Cd, Ch, Ce )615 & dU , Cd , Ch , Ce ) 586 616 !!---------------------------------------------------------------------- 587 617 !! *** ROUTINE turb_core *** … … 596 626 !! are provided at the same height 'zzu'! 597 627 !! 598 !! References : 599 !! Large & Yeager, 2004 : ??? 600 !! History : 601 !! ! XX-XX (??? ) Original code 602 !! 9.0 ! 05-08 (L. Brodeau) Rewriting and optimization 628 !! References : Large & Yeager, 2004 : ??? 603 629 !!---------------------------------------------------------------------- 604 REAL(wp), INTENT(in) :: zu ! altitude of wind measurement [m] 605 REAL(wp), INTENT(in), DIMENSION(jpi,jpj) :: & 606 sst, & ! sea surface temperature [Kelvin] 607 T_a, & ! potential air temperature [Kelvin] 608 q_sat, & ! sea surface specific humidity [kg/kg] 609 q_a, & ! specific air humidity [kg/kg] 610 dU ! wind module |U(zu)-U(0)| [m/s] 611 REAL(wp), intent(out), DIMENSION(jpi,jpj) :: & 612 Cd, & ! transfert coefficient for momentum (tau) 613 Ch, & ! transfert coefficient for temperature (Q_sens) 614 Ce ! transfert coefficient for evaporation (Q_lat) 615 616 !! * Local declarations 617 REAL(wp), DIMENSION(jpi,jpj) :: & 618 dU10, & ! dU [m/s] 619 dT, & ! air/sea temperature differeence [K] 620 dq, & ! air/sea humidity difference [K] 621 Cd_n10, & ! 10m neutral drag coefficient 622 Ce_n10, & ! 10m neutral latent coefficient 623 Ch_n10, & ! 10m neutral sensible coefficient 624 sqrt_Cd_n10, & ! root square of Cd_n10 625 sqrt_Cd, & ! root square of Cd 626 T_vpot, & ! virtual potential temperature [K] 627 T_star, & ! turbulent scale of tem. fluct. 628 q_star, & ! turbulent humidity of temp. fluct. 629 U_star, & ! turb. scale of velocity fluct. 630 L, & ! Monin-Obukov length [m] 631 zeta, & ! stability parameter at height zu 632 U_n10, & ! neutral wind velocity at 10m [m] 633 xlogt, xct, zpsi_h, zpsi_m 630 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released 631 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 coefficient 635 USE wrk_nemo, ONLY: Ce_n10 => wrk_2d_18 ! 10m neutral latent coefficient 636 USE wrk_nemo, ONLY: Ch_n10 => wrk_2d_19 ! 10m neutral sensible coefficient 637 USE wrk_nemo, ONLY: sqrt_Cd_n10 => wrk_2d_20 ! root square of Cd_n10 638 USE wrk_nemo, ONLY: sqrt_Cd => wrk_2d_21 ! root square of Cd 639 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 zu 645 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_32 648 USE wrk_nemo, ONLY: stab => iwrk_2d_1 ! 1st guess stability test integer 649 ! 650 REAL(wp) , INTENT(in ) :: zu ! altitude of wind measurement [m] 651 REAL(wp), DIMENSION(:,:), INTENT(in ) :: sst ! sea surface temperature [Kelvin] 652 REAL(wp), DIMENSION(:,:), INTENT(in ) :: T_a ! potential air temperature [Kelvin] 653 REAL(wp), DIMENSION(:,:), INTENT(in ) :: q_sat ! sea surface specific humidity [kg/kg] 654 REAL(wp), DIMENSION(:,:), INTENT(in ) :: q_a ! specific air humidity [kg/kg] 655 REAL(wp), DIMENSION(:,:), INTENT(in ) :: dU ! wind module |U(zu)-U(0)| [m/s] 656 REAL(wp), DIMENSION(:,:), INTENT( out) :: Cd ! transfert coefficient for momentum (tau) 657 REAL(wp), DIMENSION(:,:), INTENT( out) :: Ch ! transfert coefficient for temperature (Q_sens) 658 REAL(wp), DIMENSION(:,:), INTENT( out) :: Ce ! transfert coefficient for evaporation (Q_lat) 634 659 !! 635 660 INTEGER :: j_itt 636 INTEGER, PARAMETER :: nb_itt = 3 637 INTEGER, DIMENSION(jpi,jpj) :: & 638 stab ! 1st guess stability test integer 639 640 REAL(wp), PARAMETER :: & 641 grav = 9.8, & ! gravity 642 kappa = 0.4 ! von Karman s constant 661 INTEGER , PARAMETER :: nb_itt = 3 662 REAL(wp), PARAMETER :: grav = 9.8 ! gravity 663 REAL(wp), PARAMETER :: kappa = 0.4 ! von Karman s constant 643 664 !!---------------------------------------------------------------------- 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 672 644 673 !! * Start 645 674 !! Air/sea differences … … 672 701 673 702 !! Stability parameters : 674 zeta = zu/L ; zeta = sign( min(abs(zeta),10.0), zeta )675 zpsi_h = psi_h(zeta)676 zpsi_m = psi_m(zeta)703 zeta = zu/L ; zeta = sign( min(abs(zeta),10.0), zeta ) 704 zpsi_h = psi_h(zeta) 705 zpsi_m = psi_m(zeta) 677 706 678 707 !! Shifting the wind speed to 10m and neutral stability : … … 701 730 END DO 702 731 !! 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') 737 ! 703 738 END SUBROUTINE TURB_CORE_1Z 704 739 … … 717 752 !! whereas wind (dU) is at 10m. 718 753 !! 719 !! References : 720 !! Large & Yeager, 2004 : ??? 721 !! History : 722 !! 9.0 ! 06-12 (L. Brodeau) Original code for 2Z 754 !! References : Large & Yeager, 2004 : ??? 723 755 !!---------------------------------------------------------------------- 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_1 ! dU [m/s] 758 USE wrk_nemo, ONLY: dT => wrk_2d_2 ! air/sea temperature difference [K] 759 USE wrk_nemo, ONLY: dq => wrk_2d_3 ! air/sea humidity difference [K] 760 USE wrk_nemo, ONLY: Cd_n10 => wrk_2d_4 ! 10m neutral drag coefficient 761 USE wrk_nemo, ONLY: Ce_n10 => wrk_2d_5 ! 10m neutral latent coefficient 762 USE wrk_nemo, ONLY: Ch_n10 => wrk_2d_6 ! 10m neutral sensible coefficient 763 USE wrk_nemo, ONLY: sqrt_Cd_n10 => wrk_2d_7 ! root square of Cd_n10 764 USE wrk_nemo, ONLY: sqrt_Cd => wrk_2d_8 ! root square of Cd 765 USE wrk_nemo, ONLY: T_vpot => wrk_2d_9 ! virtual potential temperature [K] 766 USE wrk_nemo, ONLY: T_star => wrk_2d_10 ! turbulent scale of tem. fluct. 767 USE wrk_nemo, ONLY: q_star => wrk_2d_11 ! turbulent humidity of temp. fluct. 768 USE wrk_nemo, ONLY: U_star => wrk_2d_12 ! turb. scale of velocity fluct. 769 USE wrk_nemo, ONLY: L => wrk_2d_13 ! Monin-Obukov length [m] 770 USE wrk_nemo, ONLY: zeta_u => wrk_2d_14 ! stability parameter at height zu 771 USE wrk_nemo, ONLY: zeta_t => wrk_2d_15 ! stability parameter at height zt 772 USE wrk_nemo, ONLY: U_n10 => wrk_2d_16 ! neutral wind velocity at 10m [m] 773 USE wrk_nemo, ONLY: xlogt => wrk_2d_17, xct => wrk_2d_18, zpsi_hu => wrk_2d_19, zpsi_ht => wrk_2d_20, zpsi_m => wrk_2d_21 774 USE wrk_nemo, ONLY: stab => iwrk_2d_1 ! 1st guess stability test integer 775 !! 724 776 REAL(wp), INTENT(in) :: & 725 777 zt, & ! height for T_zt and q_zt [m] … … 738 790 q_zu ! spec. hum. shifted at zu [kg/kg] 739 791 740 !! * Local declarations741 REAL(wp), DIMENSION(jpi,jpj) :: &742 dU10, & ! dU [m/s]743 dT, & ! air/sea temperature differeence [K]744 dq, & ! air/sea humidity difference [K]745 Cd_n10, & ! 10m neutral drag coefficient746 Ce_n10, & ! 10m neutral latent coefficient747 Ch_n10, & ! 10m neutral sensible coefficient748 sqrt_Cd_n10, & ! root square of Cd_n10749 sqrt_Cd, & ! root square of Cd750 T_vpot_u, & ! virtual potential temperature [K]751 T_star, & ! turbulent scale of tem. fluct.752 q_star, & ! turbulent humidity of temp. fluct.753 U_star, & ! turb. scale of velocity fluct.754 L, & ! Monin-Obukov length [m]755 zeta_u, & ! stability parameter at height zu756 zeta_t, & ! stability parameter at height zt757 U_n10, & ! neutral wind velocity at 10m [m]758 xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m759 760 792 INTEGER :: j_itt 761 793 INTEGER, PARAMETER :: nb_itt = 3 ! number of itterations 762 INTEGER, DIMENSION(jpi,jpj) :: &763 & stab ! 1st stability test integer764 794 REAL(wp), PARAMETER :: & 765 795 grav = 9.8, & ! gravity … … 767 797 !!---------------------------------------------------------------------- 768 798 !! * Start 799 800 IF( wrk_in_use(2, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21) .OR. & 801 iwrk_in_use(2, 1) ) THEN 802 CALL ctl_stop('TURB_CORE_2Z: requested workspace arrays unavailable') ; RETURN 803 END IF 769 804 770 805 !! Initial air/sea differences … … 789 824 DO j_itt=1, nb_itt 790 825 dT = T_zu - sst ; dq = q_zu - q_sat ! Updating air/sea differences 791 T_vpot _u= T_zu*(1. + 0.608*q_zu) ! Updating virtual potential temperature at zu826 T_vpot = T_zu*(1. + 0.608*q_zu) ! Updating virtual potential temperature at zu 792 827 U_star = sqrt_Cd*dU10 ! Updating turbulent scales : (L & Y eq. (7)) 793 828 T_star = Ch/sqrt_Cd*dT ! … … 795 830 !! 796 831 L = (U_star*U_star) & ! Estimate the Monin-Obukov length at height zu 797 & / (kappa*grav/T_vpot _u*(T_star*(1.+0.608*q_zu) + 0.608*T_zu*q_star))832 & / (kappa*grav/T_vpot*(T_star*(1.+0.608*q_zu) + 0.608*T_zu*q_star)) 798 833 !! Stability parameters : 799 834 zeta_u = zu/L ; zeta_u = sign( min(abs(zeta_u),10.0), zeta_u ) … … 841 876 END DO 842 877 !! 878 IF( wrk_not_released(2, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21) .OR. & 879 iwrk_not_released(2, 1) ) CALL ctl_stop('TURB_CORE_2Z: requested workspace arrays unavailable') 880 ! 843 881 END SUBROUTINE TURB_CORE_2Z 844 882 845 883 846 884 FUNCTION psi_m(zta) !! Psis, L & Y eq. (8c), (8d), (8e) 885 !------------------------------------------------------------------------------- 886 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 887 USE wrk_nemo, ONLY: X2 => wrk_2d_33 888 USE wrk_nemo, ONLY: X => wrk_2d_34 889 USE wrk_nemo, ONLY: stabit => wrk_2d_35 890 !! 847 891 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 848 892 849 893 REAL(wp), PARAMETER :: pi = 3.141592653589793_wp 850 894 REAL(wp), DIMENSION(jpi,jpj) :: psi_m 851 REAL(wp), DIMENSION(jpi,jpj) :: X2, X, stabit 895 !------------------------------------------------------------------------------- 896 897 IF( wrk_in_use(2, 33,34,35) ) THEN 898 CALL ctl_stop('psi_m: requested workspace arrays unavailable') ; RETURN 899 ENDIF 900 852 901 X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.0) ; X = sqrt(X2) 853 902 stabit = 0.5 + sign(0.5,zta) 854 psi_m = -5.*zta*stabit & ! Stable 855 & + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable 903 psi_m = -5.*zta*stabit & ! Stable 904 & + (1. - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable 905 906 IF( wrk_not_released(2, 33,34,35) ) CALL ctl_stop('psi_m: failed to release workspace arrays') 907 ! 856 908 END FUNCTION psi_m 857 909 858 FUNCTION psi_h(zta) !! Psis, L & Y eq. (8c), (8d), (8e) 859 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 860 861 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 862 REAL(wp), DIMENSION(jpi,jpj) :: X2, X, stabit 910 911 FUNCTION psi_h( zta ) !! Psis, L & Y eq. (8c), (8d), (8e) 912 !------------------------------------------------------------------------------- 913 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 914 USE wrk_nemo, ONLY: X2 => wrk_2d_33 915 USE wrk_nemo, ONLY: X => wrk_2d_34 916 USE wrk_nemo, ONLY: stabit => wrk_2d_35 917 ! 918 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: zta 919 ! 920 REAL(wp), DIMENSION(jpi,jpj) :: psi_h 921 !------------------------------------------------------------------------------- 922 923 IF( wrk_in_use(2, 33,34,35) ) THEN 924 CALL ctl_stop('psi_h: requested workspace arrays unavailable') ; RETURN 925 ENDIF 926 863 927 X2 = sqrt(abs(1. - 16.*zta)) ; X2 = max(X2 , 1.) ; X = sqrt(X2) 864 928 stabit = 0.5 + sign(0.5,zta) 865 929 psi_h = -5.*zta*stabit & ! Stable 866 & + (1. - stabit)*(2.*log( (1. + X2)/2. )) ! Unstable 930 & + (1. - stabit)*(2.*log( (1. + X2)/2. )) ! Unstable 931 932 IF( wrk_not_released(2, 33,34,35) ) CALL ctl_stop('psi_h: failed to release workspace arrays') 933 ! 867 934 END FUNCTION psi_h 868 935
Note: See TracChangeset
for help on using the changeset viewer.
