- Timestamp:
- 2016-06-21T16:25:51+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
r6723 r6727 9 9 !! 2.0 ! 2005-04 (L. Brodeau, A.M. Treguier) additions: 10 10 !! - new bulk routine for efficiency 11 !! - WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files !!!!11 !! - WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files 12 12 !! - file names and file characteristics in namelist 13 13 !! - Implement reading of 6-hourly fields … … 18 18 !! 3.4 ! 2011-11 (C. Harris) Fill arrays required by CICE 19 19 !! 3.7 ! 2014-06 (L. Brodeau) simplification and optimization of CORE bulk 20 !! 21 !! 4.0 ! 2016-06 (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to 22 !! the CORE algorithm anymore 23 !! => based on AeroBulk (http://aerobulk.sourceforge.net/) 20 !! 4.0 ! 2016-06 (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore 21 !! ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 24 22 !!---------------------------------------------------------------------- 25 23 26 24 !!---------------------------------------------------------------------- 27 !! sbc_blk : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea) 28 !! blk_oce : computes momentum, heat and freshwater fluxes over ocean 29 !! blk_ice : computes momentum, heat and freshwater fluxes over ice 30 !! rho_air : density of (moist) air (depends on T_air, q_air and SLP 31 !! cp_air : specific heat of (moist) air (depends spec. hum. q_air) 32 !! q_sat : saturation humidity as a unction of SLP and temperature 25 !! sbc_blk : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea) 26 !! blk_oce : computes momentum, heat and freshwater fluxes over ocean 27 !! blk_ice : computes momentum, heat and freshwater fluxes over ice 28 !! rho_air : density of (moist) air (depends on T_air, q_air and SLP 29 !! cp_air : specific heat of (moist) air (depends spec. hum. q_air) 30 !! q_sat : saturation humidity as a function of SLP and temperature 31 !! L_vap : latent heat of vaporization of water as a function of temperature 33 32 !!---------------------------------------------------------------------- 34 33 USE oce ! ocean dynamics and tracers … … 49 48 USE par_ice_2 ! LIM-2 parameters 50 49 #endif 51 52 50 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 53 51 USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003) 54 52 USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013) 55 53 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31) 56 54 ! 57 55 USE iom ! I/O manager library 58 56 USE in_out_manager ! I/O manager … … 66 64 PRIVATE 67 65 68 PUBLIC sbc_blk 66 PUBLIC sbc_blk ! routine called in sbcmod module 69 67 #if defined key_lim2 || defined key_lim3 70 PUBLIC blk_ice_tau 71 PUBLIC blk_ice_flx 68 PUBLIC blk_ice_tau ! routine called in sbc_ice_lim module 69 PUBLIC blk_ice_flx ! routine called in sbc_ice_lim module 72 70 #endif 73 71 74 !!Lolo: should ultimately be moved in the module with all physical constants ? 72 !!Lolo: should ultimately be moved in the module with all physical constants ? 73 !!gm : In principle, yes. 75 74 REAL(wp), PARAMETER :: Cp_dry = 1005.0 !: Specic heat of dry air, constant pressure [J/K/kg] 76 75 REAL(wp), PARAMETER :: Cp_vap = 1860.0 !: Specic heat of water vapor, constant pressure [J/K/kg] … … 381 380 382 381 ! ... specific humidity at SST and IST tmask( 383 zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1))382 zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 384 383 !! 385 384 !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 386 385 !! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 387 386 !! (since reanalysis products provide T at z, not theta !) 388 ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1)) * rn_zqt 389 390 SELECT CASE ( nblk ) 391 CASE( np_NCAR ) ; WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008)' 392 CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003)' 393 CASE( np_COARE_3p5 ) ; WRITE(numout,*) ' "COARE 3.5" algorithm (Edson et al. 2013)' 394 CASE( np_ECMWF ) ; WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 31)' 395 END SELECT 387 ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt 396 388 397 389 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 398 399 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & 400 &Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )401 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & 402 &Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )403 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & 404 &Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )405 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & 406 &Cd, Ch, Ce, zt_zu, zq_zu, zU_zu )390 ! 391 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2 392 & Cd, Ch, Ce, zt_zu, zq_zu, zU_zu ) 393 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0 394 & Cd, Ch, Ce, zt_zu, zq_zu, zU_zu ) 395 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5 396 & Cd, Ch, Ce, zt_zu, zq_zu, zU_zu ) 397 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF 398 & Cd, Ch, Ce, zt_zu, zq_zu, zU_zu ) 407 399 CASE DEFAULT 408 400 CALL ctl_stop( 'STOP', 'sbc_blk: non-existing bulk formula selected' ) 409 401 END SELECT 410 402 411 ! Computing true air density : 412 IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN 413 !! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 414 zrhoa = rho_air(zt_zu(:,:), zq_zu(:,:), sf(jp_slp)%fnow(:,:,1)) 415 ELSE 416 !! At zt: 417 zrhoa = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 403 ! ! Compute true air density : 404 IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 405 zrhoa(:,:) = rho_air( zt_zu(:,:) , zq_zu(:,:) , sf(jp_slp)%fnow(:,:,1) ) 406 ELSE ! At zt: 407 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 418 408 END IF 419 409 420 ! ... tau module, i and j component 421 DO jj = 1, jpj 410 DO jj = 1, jpj ! tau module, i and j component 422 411 DO ji = 1, jpi 423 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd(ji,jj) ! using bulk wind speed412 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd(ji,jj) ! using bulk wind speed 424 413 taum (ji,jj) = zztmp * wndm (ji,jj) 425 414 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) … … 428 417 END DO 429 418 430 ! ... add the HF tau contribution to the wind stress module? 431 IF( lhftau ) THEN 432 taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 433 ENDIF 419 ! ! add the HF tau contribution to the wind stress module 420 IF( lhftau ) taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 421 434 422 CALL iom_put( "taum_oce", taum ) ! output wind stress module 435 423 … … 466 454 ENDIF 467 455 468 zqla(:,:) = L vap(zst(:,:)) * zevap(:,:) ! Latent Heat flux456 zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:) ! Latent Heat flux 469 457 470 458 … … 641 629 !! caution : the net upward water flux has with mm/day unit 642 630 !!--------------------------------------------------------------------- 643 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature 644 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 645 !! 646 INTEGER :: ji, jj, jl ! dummy loop indices 647 REAL(wp) :: zst2, zst3 648 REAL(wp) :: zcoef_dqlw, zcoef_dqla 649 REAL(wp) :: zztmp, z1_lsub ! temporary variable 650 !! 651 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice 652 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice 653 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqlw ! long wave heat sensitivity over ice 654 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice 655 REAL(wp), DIMENSION(:,:) , POINTER :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 631 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature 632 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 633 !! 634 INTEGER :: ji, jj, jl ! dummy loop indices 635 REAL(wp) :: zst2, zst3 ! local variable 636 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 637 REAL(wp) :: zztmp, z1_lsub ! - - 638 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice 639 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice 640 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqlw ! long wave heat sensitivity over ice 641 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice 642 REAL(wp), DIMENSION(:,:) , POINTER :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 656 643 REAL(wp), DIMENSION(:,:) , POINTER :: zrhoa 657 644 !!--------------------------------------------------------------------- … … 666 653 zcoef_dqla = -Ls * Cice * 11637800. * (-5897.8) 667 654 ! 668 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))655 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 669 656 ! 670 657 zztmp = 1. / ( 1. - albo ) … … 782 769 783 770 END SUBROUTINE blk_ice_flx 771 784 772 #endif 785 786 773 787 774 FUNCTION rho_air( ptak, pqa, pslp ) 788 775 !!------------------------------------------------------------------------------- 789 !! ** Purpose : compute density of (moist) air with eq. of state 776 !! *** FUNCTION rho_air *** 777 !! 778 !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 790 779 !! 791 780 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 792 781 !!------------------------------------------------------------------------------- 793 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak, & !: air temperature[K]794 & pqa, & !: air spec. hum.[kg/kg]795 & pslp !: pressure in[Pa]796 REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! :[kg/m^3]782 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 783 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] 784 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa] 785 REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3] 797 786 !!------------------------------------------------------------------------------- 798 787 ! 799 rho_air = pslp /(R_dry*ptak*(1._wp + rctv0*pqa))788 rho_air = pslp / ( R_dry*ptak * ( 1._wp + rctv0*pqa ) ) 800 789 ! 801 790 END FUNCTION rho_air … … 804 793 FUNCTION cp_air( pqa ) 805 794 !!------------------------------------------------------------------------------- 795 !! *** FUNCTION cp_air *** 796 !! 806 797 !! ** Purpose : Compute specific heat (Cp) of moist air 807 798 !! 808 799 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 809 800 !!------------------------------------------------------------------------------- 810 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: air spec. hum.[kg/kg]811 REAL(wp), DIMENSION(jpi,jpj) :: cp_air !:[J/K/kg]801 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] 802 REAL(wp), DIMENSION(jpi,jpj) :: cp_air ! specific heat of moist air [J/K/kg] 812 803 !!------------------------------------------------------------------------------- 813 804 ! 814 Cp_air = Cp_dry + Cp_vap *pqa805 Cp_air = Cp_dry + Cp_vap * pqa 815 806 ! 816 807 END FUNCTION cp_air … … 819 810 FUNCTION q_sat( ptak, pslp ) 820 811 !!---------------------------------------------------------------------------------- 812 !! *** FUNCTION q_sat *** 813 !! 821 814 !! ** Purpose : Specific humidity at saturation in [kg/kg] 822 815 !! Based on accurate estimate of "e_sat" … … 825 818 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 826 819 !!---------------------------------------------------------------------------------- 827 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: & 828 & ptak, & !: air temperature [K] 829 & pslp !: sea level atmospheric pressure [Pa] 830 REAL(wp), DIMENSION(jpi,jpj) :: q_sat 820 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 821 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa] 822 REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg] 831 823 ! 832 824 INTEGER :: ji, jj ! dummy loop indices 833 REAL(wp) :: ze_sat, ztmp! local scalar825 REAL(wp) :: ze_sat, ztmp ! local scalar 834 826 !!---------------------------------------------------------------------------------- 835 827 ! … … 837 829 DO ji = 1, jpi 838 830 ! 839 ztmp = rt0 /ptak(ji,jj)831 ztmp = rt0 / ptak(ji,jj) 840 832 ! 841 ! Vapour pressure at saturation [hPa] : 842 ! WMO, (Goff, 1957) 843 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) & 844 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) & 833 ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 834 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) & 835 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) & 845 836 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 ) 846 !847 q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat ) 837 ! 838 q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat ) ! 0.01 because SLP is in [Pa] 848 839 ! 849 840 END DO … … 855 846 FUNCTION gamma_moist( ptak, pqa ) 856 847 !!---------------------------------------------------------------------------------- 848 !! *** FUNCTION gamma_moist *** 849 !! 857 850 !! ** Purpose : Compute the moist adiabatic lapse-rate. 858 851 !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate … … 861 854 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 862 855 !!---------------------------------------------------------------------------------- 863 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist 864 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak, pqa ! air temperature (K) and specific humidity (kg/kg) 856 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 857 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg] 858 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate 865 859 ! 866 860 INTEGER :: ji, jj ! dummy loop indices … … 871 865 DO ji = 1, jpi 872 866 zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 873 ziRT = 1. /(R_dry*ptak(ji,jj)) ! 1/RT867 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT 874 868 gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) ) 875 869 END DO … … 879 873 880 874 881 FUNCTION L vap( psst )875 FUNCTION L_vap( psst ) 882 876 !!--------------------------------------------------------------------------------- 877 !! *** FUNCTION L_vap *** 878 !! 883 879 !! ** Purpose : Compute the latent heat of vaporization of water from temperature 884 880 !! 885 881 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 886 882 !!---------------------------------------------------------------------------------- 887 REAL(wp), DIMENSION(jpi,jpj) :: Lvap !:[J/kg]888 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst !: water temperature[K]883 REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg] 884 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K] 889 885 !!---------------------------------------------------------------------------------- 890 886 ! 891 L vap = (2.501 - 0.00237*(psst(:,:) - rt0))*1.E6892 ! 893 END FUNCTION L vap887 L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6 888 ! 889 END FUNCTION L_vap 894 890 895 891 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
