Changeset 11111 for NEMO/branches
- Timestamp:
- 2019-06-14T15:55:28+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE
- Files:
-
- 2 added
- 7 edited
-
DOM/phycst.F90 (modified) (2 diffs)
-
SBC/sbc_oce.F90 (modified) (1 diff)
-
SBC/sbcblk.F90 (modified) (29 diffs)
-
SBC/sbcblk_algo_coare.F90 (modified) (21 diffs)
-
SBC/sbcblk_algo_coare3p5.F90 (modified) (7 diffs)
-
SBC/sbcblk_algo_ecmwf.F90 (modified) (22 diffs)
-
SBC/sbcblk_algo_ncar.F90 (modified) (6 diffs)
-
SBC/sbcblk_phymbl.F90 (added)
-
SBC/sbcblk_skin.F90 (added)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/DOM/phycst.F90
r10068 r11111 50 50 REAL(wp), PUBLIC :: sice = 6.0_wp !: salinity of ice (for pisces) [psu] 51 51 REAL(wp), PUBLIC :: soce = 34.7_wp !: salinity of sea (for pisces and isf) [psu] 52 REAL(wp), PUBLIC :: rLevap = 2. 5e+6_wp !: latent heat of evaporation (water)52 REAL(wp), PUBLIC :: rLevap = 2.46e+6_wp !: latent heat of vaporization for sea-water [J/kg] #LB 53 53 REAL(wp), PUBLIC :: vkarmn = 0.4_wp !: von Karman constant 54 54 REAL(wp), PUBLIC :: stefan = 5.67e-8_wp !: Stefan-Boltzmann constant … … 66 66 REAL(wp), PUBLIC :: r1_rhos !: 1 / rhos 67 67 REAL(wp), PUBLIC :: r1_rcpi !: 1 / rcpi 68 69 70 !#LB: 71 !! Mainly used in OCE/SBC (removed from sbcblk.F90) 72 REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp !: Specic heat of dry air, constant pressure [J/K/kg] 73 REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp !: Specic heat of water vapor, constant pressure [J/K/kg] 74 REAL(wp), PARAMETER, PUBLIC :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg] 75 REAL(wp), PARAMETER, PUBLIC :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg] 76 REAL(wp), PARAMETER, PUBLIC :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622 77 REAL(wp), PARAMETER, PUBLIC :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 78 REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp !: specific heat of air (only used for ice fluxes now...) 79 REAL(wp), PARAMETER, PUBLIC :: rCd_ice = 1.4e-3_wp !: transfer coefficient over ice 80 REAL(wp), PARAMETER, PUBLIC :: albo = 0.066_wp !: ocean albedo assumed to be constant 81 ! 82 REAL(wp), PARAMETER, PUBLIC :: rho0_a = 1.2_wp !: Approx. of density of air [kg/m^3] 83 REAL(wp), PARAMETER, PUBLIC :: rho0_w = 1025._wp !: Density of sea-water (ECMWF->1025) [kg/m^3] 84 REAL(wp), PARAMETER, PUBLIC :: rCp0_w = 4190._wp !: Specific heat capacity of seawater (ECMWF 4190) [J/K/kg] 85 REAL(wp), PARAMETER, PUBLIC :: rnu0_w = 1.e-6_wp !: kinetic viscosity of water [m^2/s] 86 REAL(wp), PARAMETER, PUBLIC :: rk0_w = 0.6_wp !: thermal conductivity of water (at 20C) [W/m/K] 87 ! 88 REAL(wp), PARAMETER, PUBLIC :: emiss_w = 1._wp !: Surface emissivity (black-body long-wave radiation) of sea-water [] 89 ! !: Theoretically close to 0.97! Yet, taken equal as 1 to account for 90 ! !: the small fraction of downwelling shortwave reflected at the 91 ! !: surface (Lind & Katsaros, 1986) 92 REAL(wp), PARAMETER, PUBLIC :: rdct_qsat_salt = 0.98_wp !: reduction factor on specific humidity at saturation (q_sat(T_s)) due to salt 93 !#LB. 94 95 96 68 97 !!---------------------------------------------------------------------- 69 98 !! NEMO/OCE 4.0 , NEMO Consortium (2018) -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbc_oce.F90
r10882 r11111 149 149 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_m !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-] 150 150 151 !!---------------------------------------------------------------------- 152 !! Cool-skin/Warm-layer 153 !!---------------------------------------------------------------------- 154 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tsk !: sea-surface skin temperature out of the cool-skin/warm-layer parameterization [Celsius] 155 156 151 157 !! * Substitutions 152 158 # include "vectopt_loop_substitute.h90" -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90
r10535 r11111 15 15 !! 3.7 ! 2014-06 (L. Brodeau) simplification and optimization of CORE bulk 16 16 !! 4.0 ! 2016-06 (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore 17 !! ! ==> based on AeroBulk (http ://aerobulk.sourceforge.net/)17 !! ! ==> based on AeroBulk (https://github.com/brodeau/aerobulk/) 18 18 !! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine 19 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) … … 24 24 !! sbc_blk : bulk formulation as ocean surface boundary condition 25 25 !! blk_oce : computes momentum, heat and freshwater fluxes over ocean 26 !! rho_air : density of (moist) air (depends on T_air, q_air and SLP27 !! cp_air : specific heat of (moist) air (depends spec. hum. q_air)28 !! q_sat : saturation humidity as a function of SLP and temperature29 !! L_vap : latent heat of vaporization of water as a function of temperature30 26 !! sea-ice case only : 31 27 !! blk_ice_tau : provide the air-ice stress … … 60 56 USE prtctl ! Print control 61 57 58 USE sbcblk_phymbl !LB: all thermodynamics functions, rho_air, q_sat, etc... 59 60 62 61 IMPLICIT NONE 63 62 PRIVATE … … 70 69 PUBLIC blk_ice_qcn ! routine called in icesbc 71 70 #endif 72 73 !!Lolo: should ultimately be moved in the module with all physical constants ?74 !!gm : In principle, yes.75 REAL(wp), PARAMETER :: Cp_dry = 1005.0 !: Specic heat of dry air, constant pressure [J/K/kg]76 REAL(wp), PARAMETER :: Cp_vap = 1860.0 !: Specic heat of water vapor, constant pressure [J/K/kg]77 REAL(wp), PARAMETER :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg]78 REAL(wp), PARAMETER :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg]79 REAL(wp), PARAMETER :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.62280 REAL(wp), PARAMETER :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.60881 71 82 72 INTEGER , PARAMETER :: jpfld =10 ! maximum number of files to read … … 94 84 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read) 95 85 96 ! !!! Bulk parameters97 REAL(wp), PARAMETER :: cpa = 1000.5 ! specific heat of air (only used for ice fluxes now...)98 REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation99 REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant100 REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! transfer coefficient over ice101 REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant102 !103 86 ! !!* Namelist namsbc_blk : bulk parameters 104 87 LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008) … … 124 107 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 125 108 109 !LB: 110 LOGICAL :: ln_skin ! use the cool-skin / warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 111 !LB: => sbc_oce.F90: REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: tsk ! sea-surface skin temperature out of the cool-skin/warm-layer parameterization [Celsius] 112 !LB. 113 126 114 INTEGER :: nblk ! choice of the bulk algorithm 127 115 ! ! associated indices: … … 150 138 IF( sbc_blk_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' ) 151 139 END FUNCTION sbc_blk_alloc 140 141 !LB: 142 INTEGER FUNCTION sbc_blk_cswl_alloc() 143 !!------------------------------------------------------------------- 144 !! *** ROUTINE sbc_blk_cswl_alloc *** 145 !!------------------------------------------------------------------- 146 !PRINT *, '*** LB: allocating tsk!' 147 ALLOCATE( tsk(jpi,jpj), STAT=sbc_blk_cswl_alloc ) 148 !PRINT *, '*** LB: done!' 149 CALL mpp_sum ( 'sbcblk', sbc_blk_cswl_alloc ) 150 IF( sbc_blk_cswl_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_cswl_alloc: failed to allocate arrays' ) 151 END FUNCTION sbc_blk_cswl_alloc 152 !LB. 152 153 153 154 … … 173 174 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm 174 175 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 175 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 176 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, & 177 & ln_skin ! cool-skin / warm-layer param. !LB 176 178 !!--------------------------------------------------------------------- 177 179 ! … … 198 200 IF( ln_ECMWF ) THEN ; nblk = np_ECMWF ; ioptio = ioptio + 1 ; ENDIF 199 201 ! 202 !LB: 203 ! !** initialization of the cool-skin / warm-layer parametrization 204 IF( ln_skin ) THEN 205 IF ( ln_NCAR ) CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' ) 206 ! ! allocate array(s) for cool-skin/warm-layer param. 207 IF( sbc_blk_cswl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 208 END IF 209 !LB. 210 200 211 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 201 212 ! … … 279 290 END SELECT 280 291 ! 292 WRITE(numout,*) 293 WRITE(numout,*) ' use cool-skin/warm-layer parameterization (SSST) ln_skin = ', ln_skin !LB 294 ! 281 295 ENDIF 282 296 ! … … 413 427 414 428 ! ----------------------------------------------------------------------------- ! 415 ! I Radiative FLUXES!429 ! I Solar FLUX ! 416 430 ! ----------------------------------------------------------------------------- ! 417 431 … … 422 436 ENDIF 423 437 424 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave425 438 426 439 ! ----------------------------------------------------------------------------- ! … … 429 442 430 443 ! ... specific humidity at SST and IST tmask( 431 zsq(:,:) = 0.98* q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )444 zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 432 445 !! 433 446 !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate … … 444 457 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5 445 458 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 446 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF 447 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 459 460 !LB: Skin!!! 461 CASE( np_ECMWF ) 462 IF ( ln_skin ) THEN 463 !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL optional arrays!!!' 464 !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => BEFORE ZST(40:50,30) =', zst(40:50,30) 465 CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF 466 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 467 & Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1)) 468 !LB: "zst" and "zsq" have been updated with skin temp. !!! (from bulk SST)... 469 zst(:,:) = zst(:,:)*tmask(:,:,1) 470 zsq(:,:) = zsq(:,:)*tmask(:,:,1) 471 !IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => AFTER ZST(40:50,30) =', zst(40:50,30) 472 ELSE 473 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!' 474 CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF 475 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 476 END IF 477 !LB. 478 448 479 CASE DEFAULT 449 480 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 450 481 END SELECT 482 483 484 !LB: cool-skin/warm-layer: 485 IF ( ln_skin ) tsk(:,:) = zst(:,:) 486 451 487 452 488 ! ! Compute true air density : … … 507 543 508 544 545 ! ----------------------------------------------------------------------------- ! 546 ! III Net longwave radiative FLUX ! 547 ! ----------------------------------------------------------------------------- ! 548 549 !! LB: now after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin==.TRUE.) 550 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - emiss_w * stefan * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 551 552 509 553 IF(ln_ctl) THEN 510 554 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce : ' ) … … 519 563 520 564 ! ----------------------------------------------------------------------------- ! 521 ! I IITotal FLUXES !565 ! IV Total FLUXES ! 522 566 ! ----------------------------------------------------------------------------- ! 523 567 ! … … 527 571 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 528 572 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 529 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST 573 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST !LB??? pst is Celsius !? 530 574 & + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac & ! add liquid precip heat content at Tair 531 575 & * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & … … 539 583 #endif 540 584 ! 541 IF ( nn_ice == 0 ) THEN 585 !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN 586 CALL iom_put( "rho_air" , zrhoa ) ! output air density (kg/m^3) !#LB 542 587 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 543 588 CALL iom_put( "qsb_oce" , - zqsb ) ! output downward sensible heat over the ocean … … 551 596 CALL iom_put( 'snowpre', sprecip ) ! Snow 552 597 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 553 ENDIF 598 IF ( ln_skin ) THEN 599 CALL iom_put( "t_skin" , (zst - rt0) * tmask(:,:,1) ) ! T_skin in Celsius 600 CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) ) ! T_skin - SST temperature difference... 601 END IF 602 !!#LB. ENDIF 554 603 ! 555 604 IF(ln_ctl) THEN … … 564 613 565 614 566 567 FUNCTION rho_air( ptak, pqa, pslp )568 !!-------------------------------------------------------------------------------569 !! *** FUNCTION rho_air ***570 !!571 !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere572 !!573 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)574 !!-------------------------------------------------------------------------------575 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]576 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]577 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa]578 REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3]579 !!-------------------------------------------------------------------------------580 !581 rho_air = pslp / ( R_dry*ptak * ( 1._wp + rctv0*pqa ) )582 !583 END FUNCTION rho_air584 585 586 FUNCTION cp_air( pqa )587 !!-------------------------------------------------------------------------------588 !! *** FUNCTION cp_air ***589 !!590 !! ** Purpose : Compute specific heat (Cp) of moist air591 !!592 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)593 !!-------------------------------------------------------------------------------594 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]595 REAL(wp), DIMENSION(jpi,jpj) :: cp_air ! specific heat of moist air [J/K/kg]596 !!-------------------------------------------------------------------------------597 !598 Cp_air = Cp_dry + Cp_vap * pqa599 !600 END FUNCTION cp_air601 602 603 FUNCTION q_sat( ptak, pslp )604 !!----------------------------------------------------------------------------------605 !! *** FUNCTION q_sat ***606 !!607 !! ** Purpose : Specific humidity at saturation in [kg/kg]608 !! Based on accurate estimate of "e_sat"609 !! aka saturation water vapor (Goff, 1957)610 !!611 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)612 !!----------------------------------------------------------------------------------613 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]614 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa]615 REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg]616 !617 INTEGER :: ji, jj ! dummy loop indices618 REAL(wp) :: ze_sat, ztmp ! local scalar619 !!----------------------------------------------------------------------------------620 !621 DO jj = 1, jpj622 DO ji = 1, jpi623 !624 ztmp = rt0 / ptak(ji,jj)625 !626 ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)627 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) &628 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) &629 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 )630 !631 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]632 !633 END DO634 END DO635 !636 END FUNCTION q_sat637 638 639 FUNCTION gamma_moist( ptak, pqa )640 !!----------------------------------------------------------------------------------641 !! *** FUNCTION gamma_moist ***642 !!643 !! ** Purpose : Compute the moist adiabatic lapse-rate.644 !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate645 !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html646 !!647 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)648 !!----------------------------------------------------------------------------------649 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]650 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg]651 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate652 !653 INTEGER :: ji, jj ! dummy loop indices654 REAL(wp) :: zrv, ziRT ! local scalar655 !!----------------------------------------------------------------------------------656 !657 DO jj = 1, jpj658 DO ji = 1, jpi659 zrv = pqa(ji,jj) / (1. - pqa(ji,jj))660 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT661 gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) )662 END DO663 END DO664 !665 END FUNCTION gamma_moist666 667 668 FUNCTION L_vap( psst )669 !!---------------------------------------------------------------------------------670 !! *** FUNCTION L_vap ***671 !!672 !! ** Purpose : Compute the latent heat of vaporization of water from temperature673 !!674 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)675 !!----------------------------------------------------------------------------------676 REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg]677 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]678 !!----------------------------------------------------------------------------------679 !680 L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6681 !682 END FUNCTION L_vap683 615 684 616 #if defined key_si3 … … 710 642 ! 711 643 ! set transfer coefficients to default sea-ice values 712 Cd_atm(:,:) = Cd_ice713 Ch_atm(:,:) = Cd_ice714 Ce_atm(:,:) = Cd_ice644 Cd_atm(:,:) = rCd_ice 645 Ch_atm(:,:) = rCd_ice 646 Ce_atm(:,:) = rCd_ice 715 647 716 648 wndm_ice(:,:) = 0._wp !!gm brutal.... … … 737 669 ENDIF 738 670 739 !! CALL iom_put( " Cd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef.671 !! CALL iom_put( "rCd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef. 740 672 !! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef. 741 673 … … 792 724 REAL(wp) :: zst3 ! local variable 793 725 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 794 REAL(wp) :: zztmp, z1_rLsub 726 REAL(wp) :: zztmp, z1_rLsub ! - - 795 727 REAL(wp) :: zfr1, zfr2 ! local variables 796 728 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature … … 803 735 !!--------------------------------------------------------------------- 804 736 ! 805 zcoef_dqlw = 4. 0 * 0.95 * Stef! local scalars806 zcoef_dqla = - Ls * 11637800. * (-5897.8)737 zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars 738 zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 807 739 ! 808 740 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) … … 824 756 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 825 757 ! Long Wave (lw) 826 z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef* ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)758 z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 827 759 ! lw sensitivity 828 760 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 … … 834 766 ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 835 767 ! Sensible Heat 836 z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa* Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))768 z_qsb(ji,jj,jl) = zrhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1)) 837 769 ! Latent Heat 838 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls* Ch_atm(ji,jj) * wndm_ice(ji,jj) * &770 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * rLsub * Ch_atm(ji,jj) * wndm_ice(ji,jj) * & 839 771 & ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 840 772 ! Latent heat sensitivity for ice (Dqla/Dt) … … 847 779 848 780 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 849 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa* Ch_atm(ji,jj) * wndm_ice(ji,jj)781 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj) 850 782 851 783 ! ----------------------------! … … 1064 996 ! generic drag over a cell partly covered by ice 1065 997 !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag 1066 !! & Cd_ice * at_i_b(:,:) + & ! pure ice drag998 !! & rCd_ice * at_i_b(:,:) + & ! pure ice drag 1067 999 !! & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu ! change due to sea-ice morphology 1068 1000 1069 1001 ! ice-atm drag 1070 Cd(:,:) = Cd_ice + & ! pure ice drag1002 Cd(:,:) = rCd_ice + & ! pure ice drag 1071 1003 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology 1072 1004 … … 1141 1073 ! Atmospheric and Surface Variables 1142 1074 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin 1143 zqo_sat(:,:) = 0.98_wp* q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg]1144 zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg]1075 zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg] 1076 zqi_sat(:,:) = q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg] !LB: no 0.98 !!(rdct_qsat_salt) 1145 1077 ! 1146 1078 DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare.F90
r10069 r11111 1 1 MODULE sbcblk_algo_coare 2 2 !!====================================================================== 3 !! *** MODULE sbcblk_algo_coare *** 4 !! Computes turbulent components of surface fluxes 5 !! according to Fairall et al. 2003 (COARE v3) 6 !! 3 !! *** MODULE sbcblk_algo_coare *** 4 !! Computes: 7 5 !! * bulk transfer coefficients C_D, C_E and C_H 8 6 !! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed … … 10 8 !! => all these are used in bulk formulas in sbcblk.F90 11 9 !! 12 !! Using the bulk formulation/param. of COARE v3, Fairall et al. 2003 13 !! 10 !! Using the bulk formulation/param. of COARE v3, Fairall et al., 2003 14 11 !! 15 12 !! Routine turb_coare maintained and developed in AeroBulk 16 !! (http ://aerobulk.sourceforge.net/)13 !! (https://github.com/brodeau/aerobulk) 17 14 !! 18 !! Author: Laurent Brodeau, 2016, brodeau@gmail.com 19 !! 20 !!====================================================================== 21 !! History : 3.6 ! 2016-02 (L.Brodeau) Original code 15 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 16 !!---------------------------------------------------------------------- 17 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code 22 18 !!---------------------------------------------------------------------- 23 19 … … 27 23 !! returns the effective bulk wind speed at 10m 28 24 !!---------------------------------------------------------------------- 29 USE oce ! ocean dynamics and tracers30 USE dom_oce ! ocean space and time domain31 USE phycst ! physical constants32 USE sbc_oce ! Surface boundary condition: ocean fields33 USE sbcwave, ONLY : cdn_wave ! wave module25 USE oce ! ocean dynamics and tracers 26 USE dom_oce ! ocean space and time domain 27 USE phycst ! physical constants 28 USE sbc_oce ! Surface boundary condition: ocean fields 29 USE sbcwave, ONLY : cdn_wave ! wave module 34 30 #if defined key_si3 || defined key_cice 35 USE sbc_ice ! Surface boundary condition: ice fields31 USE sbc_ice ! Surface boundary condition: ice fields 36 32 #endif 37 !33 USE sbcblk_phymbl !LB: all thermodynamics functions, rho_air, q_sat, etc... 38 34 USE in_out_manager ! I/O manager 39 35 USE iom ! I/O manager library … … 50 46 REAL(wp), PARAMETER :: zi0 = 600._wp ! scale height of the atmospheric boundary layer... 51 47 REAL(wp), PARAMETER :: Beta0 = 1.250_wp ! gustiness parameter 52 REAL(wp), PARAMETER :: rctv0 = 0.608_wp ! constant to obtain virtual temperature... 48 49 INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations 53 50 54 51 !!---------------------------------------------------------------------- … … 61 58 !! *** ROUTINE turb_coare *** 62 59 !! 63 !! 2015: L. Brodeau (brodeau@gmail.com)64 !!65 60 !! ** Purpose : Computes turbulent transfert coefficients of surface 66 61 !! fluxes according to Fairall et al. (2003) 67 62 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 68 !! 69 !! ** Method : Monin Obukhov Similarity Theory 70 !!---------------------------------------------------------------------- 63 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 64 !! 71 65 !! 72 66 !! INPUT : … … 88 82 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 89 83 !! * q_zu : specific humidity of air // [kg/kg] 90 !! * U_blk : bulk wind at 10m [m/s] 91 !!---------------------------------------------------------------------- 84 !! * U_blk : bulk wind speed at 10m [m/s] 85 !! 86 !! 87 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 88 !!---------------------------------------------------------------------------------- 92 89 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 93 90 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] … … 106 103 ! 107 104 INTEGER :: j_itt 108 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 109 INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations 105 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 110 106 111 107 REAL(wp), DIMENSION(jpi,jpj) :: & … … 125 121 126 122 !! First guess of temperature and humidity at height zu: 127 t_zu = MAX( t_zt , 0.0)! who knows what's given on masked-continental regions...128 q_zu = MAX( q_zt , 1.e-6)! "123 t_zu = MAX( t_zt , 199.0_wp ) ! who knows what's given on masked-continental regions... 124 q_zu = MAX( q_zt , 1.e-6_wp ) ! " 129 125 130 126 !! Pot. temp. difference (and we don't want it to be 0!) 131 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6 ), dt_zu )132 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9 ), dq_zu )133 134 znu_a = visc_air(t_z t) ! Air viscosity (m^2/s) at zt given from temperature in (K)127 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 128 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 129 130 znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K) 135 131 136 132 ztmp2 = 0.5*0.5 ! initial guess for wind gustiness contribution … … 144 140 145 141 z0 = alfa_charn(U_blk)*u_star*u_star/grav + 0.11*znu_a/u_star 146 z0t = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ! WARNING: 1/z0t ! 142 z0 = MIN(ABS(z0), 0.001) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 143 z0t = 1. / ( 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 144 z0t = MIN(ABS(z0t), 0.001) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 147 145 148 146 ztmp2 = vkarmn/ztmp0 149 147 Cd = ztmp2*ztmp2 ! first guess of Cd 150 148 151 ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd 152 153 !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth 154 ztmp2 = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk) !! Ribu Bulk Richardson number 155 ztmp1 = 0.5 + sign(0.5 , ztmp2) 149 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 150 151 ztmp2 = grav*zu*(dt_zu + rctv0*t_zu*dq_zu)/(t_zu*U_blk*U_blk) !! Ribu Bulk Richardson number ; !Ribcu = -zu/(zi0*0.004*Beta0**3) !! Saturation Rib, zi0 = tropicalbound. layer depth 152 153 !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2): 154 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 156 155 ztmp0 = ztmp0*ztmp2 157 !! Ribu < 0 Ribu > 0 Beta = 1.25158 zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) &159 & + ztmp1 * (ztmp0*(1. + 27./9.*ztmp2/ztmp0))156 zeta_u = (1.-ztmp1) * (ztmp0/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & ! Ribu < 0 157 & + ztmp1 * (ztmp0*(1. + 27./9.*ztmp2/ztmp0)) ! Ribu > 0 158 !#LOLO: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 160 159 161 160 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 162 ztmp0 = vkarmn/(LOG(zu*z0t) - psi_h_coare(zeta_u))161 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u)) 163 162 164 163 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) … … 168 167 ! What's need to be done if zt /= zu: 169 168 IF( .NOT. l_zt_equal_zu ) THEN 170 169 !! First update of values at zu (or zt for wind) 171 170 zeta_t = zt*zeta_u/zu 172 173 !! First update of values at zu (or zt for wind)174 171 ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t) 175 ztmp1 = log(zt/zu) + ztmp0172 ztmp1 = LOG(zt/zu) + ztmp0 176 173 t_zu = t_zt - t_star/vkarmn*ztmp1 177 174 q_zu = q_zt - q_star/vkarmn*ztmp1 178 q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity : 179 180 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 181 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 182 175 q_zu = (0.5 + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 176 ! 177 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 178 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 183 179 END IF 184 180 … … 188 184 !!Inverse of Monin-Obukov length (1/L) : 189 185 ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star) ! 1/L == 1/[Monin-Obukhov length] 186 ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO 190 187 191 188 ztmp1 = u_star*u_star ! u*^2 … … 193 190 !! Update wind at 10m taking into acount convection-related wind gustiness: 194 191 ! Ug = Beta*w* (Beta = 1.25, Fairall et al. 2003, Eq.8): 195 ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0. ))**(2./3.) ! => ztmp2 == Ug^2192 ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2./3.) ! => ztmp2 == Ug^2 196 193 !! ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before 600. 197 U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2 ) ! include gustiness in bulk wind speed194 U_blk = MAX(sqrt(U_zu*U_zu + ztmp2), 0.2_wp) ! include gustiness in bulk wind speed 198 195 ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 199 196 … … 203 200 !! Roughness lengthes z0, z0t (z0q = z0t) : 204 201 z0 = ztmp2*ztmp1/grav + 0.11*znu_a/u_star ! Roughness length (eq.6) 205 ztmp1 = z0*u_star/znu_a 206 z0t = min( 1.1E-4 , 5.5E-5*ztmp1**(-0.6) )! Scalar roughness for both theta and q (eq.28)202 ztmp1 = z0*u_star/znu_a ! Re_r: roughness Reynolds number 203 z0t = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1**(-0.6_wp) ) ! Scalar roughness for both theta and q (eq.28) 207 204 208 205 !! Stability parameters: 209 zeta_u = zu*ztmp0 ; zeta_u = sign( min(abs(zeta_u),50.0), zeta_u ) 206 zeta_u = zu*ztmp0 207 zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u ) 210 208 IF( .NOT. l_zt_equal_zu ) THEN 211 zeta_t = zt*ztmp0 ; zeta_t = sign( min(abs(zeta_t),50.0), zeta_t ) 209 zeta_t = zt*ztmp0 210 zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t ) 212 211 END IF 213 212 214 213 !! Turbulent scales at zu=10m : 215 214 ztmp0 = psi_h_coare(zeta_u) 216 ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0)215 ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 217 216 218 217 t_star = dt_zu*ztmp1 219 218 q_star = dq_zu*ztmp1 220 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))219 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) 221 220 222 221 IF( .NOT. l_zt_equal_zu ) THEN … … 259 258 !! Wind greater than 18 m/s : alfa = 0.018 260 259 !! 261 !! Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)260 !! Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 262 261 !!------------------------------------------------------------------- 263 262 REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn … … 274 273 ! 275 274 ! Charnock's constant, increases with the wind : 276 zgt10 = 0.5 + SIGN(0.5 ,(zw - 10.)) ! If zw<10. --> 0, else --> 1277 zgt18 = 0.5 + SIGN(0.5 ,(zw - 18.)) ! If zw<18. --> 0, else --> 1275 zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1 276 zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1 278 277 ! 279 278 alfa_charn(ji,jj) = (1. - zgt10)*0.011 & ! wind is lower than 10 m/s … … 294 293 !! 295 294 !! Author: L. Brodeau, june 2016 / AeroBulk 296 !! (https:// sourceforge.net/p/aerobulk)295 !! (https://github.com/brodeau/aerobulk/) 297 296 !!------------------------------------------------------------------------ 298 297 REAL(wp), DIMENSION(jpi,jpj) :: One_on_L !: 1./(Monin Obukhov length) [m^-1] … … 330 329 !! form from Beljaars and Holtslag (1991) 331 330 !! 332 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)331 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 333 332 !!---------------------------------------------------------------------------------- 334 333 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare … … 356 355 zf = zta*zta 357 356 zf = zf/(1. + zf) 358 zc = MIN(50. , 0.35*zta)359 zstab = 0.5 + SIGN(0.5 , zta)357 zc = MIN(50._wp, 0.35_wp*zta) 358 zstab = 0.5 + SIGN(0.5_wp, zta) 360 359 ! 361 360 psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) … … 383 382 !! 384 383 !! Author: L. Brodeau, june 2016 / AeroBulk 385 !! (https:// sourceforge.net/p/aerobulk)384 !! (https://github.com/brodeau/aerobulk/) 386 385 !!---------------------------------------------------------------- 387 386 !! … … 408 407 zf = zta*zta 409 408 zf = zf/(1. + zf) 410 zc = MIN(50. ,0.35*zta)411 zstab = 0.5 + SIGN(0.5 , zta)409 zc = MIN(50._wp,0.35_wp*zta) 410 zstab = 0.5 + SIGN(0.5_wp, zta) 412 411 ! 413 412 psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & … … 420 419 END FUNCTION psi_h_coare 421 420 422 423 FUNCTION visc_air( ptak ) 424 !!--------------------------------------------------------------------- 425 !! Air kinetic viscosity (m^2/s) given from temperature in degrees... 426 !! 427 !! Author: L. Brodeau, june 2016 / AeroBulk 428 !! (https://sourceforge.net/p/aerobulk) 429 !!--------------------------------------------------------------------- 430 !! 431 REAL(wp), DIMENSION(jpi,jpj) :: visc_air 432 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K) 433 ! 434 INTEGER :: ji, jj ! dummy loop indices 435 REAL(wp) :: ztc, ztc2 ! local scalar 436 ! 437 DO jj = 1, jpj 438 DO ji = 1, jpi 439 ztc = ptak(ji,jj) - rt0 ! air temp, in deg. C 440 ztc2 = ztc*ztc 441 visc_air(ji,jj) = 1.326E-5*(1. + 6.542E-3*ztc + 8.301E-6*ztc2 - 4.84E-9*ztc2*ztc) 442 END DO 443 END DO 444 ! 445 END FUNCTION visc_air 446 447 421 !!====================================================================== 448 422 END MODULE sbcblk_algo_coare -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p5.F90
r10069 r11111 14 14 !! 15 15 !! Routine turb_coare3p5 maintained and developed in AeroBulk 16 !! (http ://aerobulk.sourceforge.net/)16 !! (https://github.com/brodeau/aerobulk/) 17 17 !! 18 18 !! Author: Laurent Brodeau, 2016, brodeau@gmail.com … … 51 51 REAL(wp), PARAMETER :: zi0 = 600. ! scale height of the atmospheric boundary layer...1 52 52 REAL(wp), PARAMETER :: Beta0 = 1.25 ! gustiness parameter 53 REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature...54 53 55 54 !!---------------------------------------------------------------------- … … 68 67 !! ** Method : Monin Obukhov Similarity Theory 69 68 !! 70 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)69 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 71 70 !! 72 71 !! INPUT : … … 265 264 !! 266 265 !! Author: L. Brodeau, june 2016 / AeroBulk 267 !! (https:// sourceforge.net/p/aerobulk)266 !! (https://github.com/brodeau/aerobulk/) 268 267 !!------------------------------------------------------------------------ 269 268 REAL(wp), DIMENSION(jpi,jpj) :: One_on_L !: 1./(Monin Obukhov length) [m^-1] … … 300 299 !! form from Beljaars and Holtslag (1991) 301 300 !! 302 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)301 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 303 302 !!---------------------------------------------------------------------------------- 304 303 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare … … 353 352 !! 354 353 !! Author: L. Brodeau, june 2016 / AeroBulk 355 !! (https:// sourceforge.net/p/aerobulk)354 !! (https://github.com/brodeau/aerobulk/) 356 355 !!---------------------------------------------------------------- 357 356 !! … … 396 395 !! 397 396 !! Author: L. Brodeau, june 2016 / AeroBulk 398 !! (https:// sourceforge.net/p/aerobulk)397 !! (https://github.com/brodeau/aerobulk/) 399 398 !!--------------------------------------------------------------------- 400 399 REAL(wp), DIMENSION(jpi,jpj) :: visc_air -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r10069 r11111 1 1 MODULE sbcblk_algo_ecmwf 2 2 !!====================================================================== 3 !! *** MODULE sbcblk_algo_ecmwf *** 4 !! Computes turbulent components of surface fluxes 5 !! according to the method in IFS of the ECMWF model 6 !! 3 !! *** MODULE sbcblk_algo_ecmwf *** 4 !! Computes: 7 5 !! * bulk transfer coefficients C_D, C_E and C_H 8 6 !! * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed … … 10 8 !! => all these are used in bulk formulas in sbcblk.F90 11 9 !! 12 !! Using the bulk formulation/param. of IFS of ECMWF (cycle 31r2)10 !! Using the bulk formulation/param. of IFS of ECMWF (cycle 40r1) 13 11 !! based on IFS doc (avaible online on the ECMWF's website) 14 12 !! 13 !! Routine turb_ecmwf maintained and developed in AeroBulk 14 !! (https://github.com/brodeau/aerobulk) 15 15 !! 16 !! Routine turb_ecmwf maintained and developed in AeroBulk 17 !! (http://aerobulk.sourceforge.net/) 18 !! 19 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 16 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk) 20 17 !!---------------------------------------------------------------------- 21 18 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code … … 41 38 42 39 USE sbc_oce ! Surface boundary condition: ocean fields 40 USE sbcblk_phymbl !LB: all thermodynamics functions, rho_air, q_sat, etc... 41 USE sbcblk_skin ! cool-skin/warm layer scheme (CSWL_ECMWF) 43 42 44 43 IMPLICIT NONE … … 52 51 REAL(wp), PARAMETER :: zi0 = 1000. ! scale height of the atmospheric boundary layer...1 53 52 REAL(wp), PARAMETER :: Beta0 = 1. ! gustiness parameter ( = 1.25 in COAREv3) 54 REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature...55 REAL(wp), PARAMETER :: Cp_dry = 1005.0 ! Specic heat of dry air, constant pressure [J/K/kg]56 REAL(wp), PARAMETER :: Cp_vap = 1860.0 ! Specic heat of water vapor, constant pressure [J/K/kg]57 53 REAL(wp), PARAMETER :: alpha_M = 0.11 ! For roughness length (smooth surface term) 58 54 REAL(wp), PARAMETER :: alpha_H = 0.40 ! (Chapter 3, p.34, IFS doc Cy31r1) 59 55 REAL(wp), PARAMETER :: alpha_Q = 0.62 ! 56 57 INTEGER , PARAMETER :: nb_itt = 5 ! number of itterations 58 59 !! Cool-Skin / Warm-Layer related parameters: 60 REAL(wp), PARAMETER :: & 61 & rdt0 = 3600.*1.5, & !: time step 62 & rd0 = 3. , & !: Depth scale [m], "d" in Eq.11 (Zeng & Beljaars 2005) 63 & rNu0 = 0.5 !: Nu (exponent of temperature profile) Eq.11 64 ! !: (Zeng & Beljaars 2005) !: set to 0.5 instead of 65 ! !: 0.3 to respect a warming of +3 K in calm 66 ! !: condition for the insolation peak of +1000W/m^2 67 INTEGER, PARAMETER :: & 68 & nb_itt_wl = 10 !: number of sub-itterations for solving the differential equation in warm-layer part 69 ! !: => use "nb_itt_wl = 1" for No itteration! => way cheaper !!! 70 ! !: => assumes balance between the last 2 terms of Eq.11 (lhs of eq.11 = 0) 71 ! !: => in that case no need for sub-itterations ! 72 ! !: => ACCEPTABLE IN MOST CONDITIONS ! (UNLESS: sunny + very calm/low-wind conditions) 73 ! !: => Otherwize use "nb_itt_wl = 10" 74 75 76 77 60 78 !!---------------------------------------------------------------------- 61 79 CONTAINS 62 80 63 SUBROUTINE TURB_ECMWF( zt, zu, sst, t_zt, ssq , q_zt , U_zu, & 64 & Cd, Ch, Ce , t_zu, q_zu, U_blk, & 65 & Cdn, Chn, Cen ) 81 SUBROUTINE TURB_ECMWF( zt, zu, T_s, t_zt, q_s, q_zt, U_zu, & 82 & Cd, Ch, Ce, t_zu, q_zu, U_blk, & 83 & Cdn, Chn, Cen, & 84 & Qsw, rad_lw, slp ) 66 85 !!---------------------------------------------------------------------------------- 67 86 !! *** ROUTINE turb_ecmwf *** 68 87 !! 69 !! 2015: L. Brodeau (brodeau@gmail.com)70 !!71 88 !! ** Purpose : Computes turbulent transfert coefficients of surface 72 !! fluxes according to IFS doc. (cycle 31)89 !! fluxes according to IFS doc. (cycle 40) 73 90 !! If relevant (zt /= zu), adjust temperature and humidity from height zt to zu 74 !! 75 !! ** Method : Monin Obukhov Similarity Theory 91 !! Returns the effective bulk wind speed at 10m to be used in the bulk formulas 92 !! 93 !! Applies the cool-skin warm-layer correction of the SST to T_s 94 !! if the net shortwave flux at the surface (Qsw), the downwelling longwave 95 !! radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp) 96 !! are provided as (optional) arguments! 76 97 !! 77 98 !! INPUT : … … 80 101 !! * zu : height for wind speed (generally 10m) [m] 81 102 !! * U_zu : scalar wind speed at 10m [m/s] 82 !! * sst : SST [K]83 103 !! * t_zt : potential air temperature at zt [K] 84 !! * ssq : specific humidity at saturation at SST [kg/kg]85 104 !! * q_zt : specific humidity of air at zt [kg/kg] 86 105 !! 106 !! INPUT/OUTPUT: 107 !! ------------- 108 !! * T_s : SST or skin temperature [K] 109 !! * q_s : SSQ aka saturation specific humidity at temp. T_s [kg/kg] 110 !! -> doesn't need to be given a value if skin temp computed (in case l_use_skin=True) 111 !! -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=False) 112 !! 113 !! OPTIONAL INPUT (will trigger l_use_skin=TRUE if present!): 114 !! --------------- 115 !! * Qsw : net solar flux (after albedo) at the surface (>0) [W/m^2] 116 !! * rad_lw : downwelling longwave radiation at the surface (>0) [W/m^2] 117 !! * slp : sea-level pressure [Pa] 87 118 !! 88 119 !! OUTPUT : … … 93 124 !! * t_zu : pot. air temperature adjusted at wind height zu [K] 94 125 !! * q_zu : specific humidity of air // [kg/kg] 95 !! * U_blk : bulk wind at 10m[m/s]96 !! 97 !! 98 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)126 !! * U_blk : bulk wind speed at 10m [m/s] 127 !! 128 !! 129 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 99 130 !!---------------------------------------------------------------------------------- 100 131 REAL(wp), INTENT(in ) :: zt ! height for t_zt and q_zt [m] 101 132 REAL(wp), INTENT(in ) :: zu ! height for U_zu [m] 102 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: sst! sea surface temperature [Kelvin]133 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: T_s ! sea surface temperature [Kelvin] 103 134 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt ! potential air temperature [Kelvin] 104 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: ssq! sea surface specific humidity [kg/kg]105 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity 135 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: q_s ! sea surface specific humidity [kg/kg] 136 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt ! specific air humidity at zt [kg/kg] 106 137 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu ! relative wind module at zu [m/s] 107 138 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cd ! transfer coefficient for momentum (tau) … … 113 144 REAL(wp), INTENT( out), DIMENSION(jpi,jpj) :: Cdn, Chn, Cen ! neutral transfer coefficients 114 145 ! 146 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: Qsw ! [W/m^2] 147 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: rad_lw ! [W/m^2] 148 REAL(wp), INTENT(in ), OPTIONAL, DIMENSION(jpi,jpj) :: slp ! [Pa] 149 ! 115 150 INTEGER :: j_itt 116 LOGICAL :: 117 INTEGER , PARAMETER :: nb_itt = 4 ! number of itterations118 !119 REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star,&151 LOGICAL :: l_zt_equal_zu = .FALSE. ! if q and t are given at same height as U 152 ! 153 REAL(wp), DIMENSION(jpi,jpj) :: & 154 & u_star, t_star, q_star, & 120 155 & dt_zu, dq_zu, & 121 156 & znu_a, & !: Nu_air, Viscosity of air 122 157 & Linv, & !: 1/L (inverse of Monin Obukhov length... 123 158 & z0, z0t, z0q 159 ! 160 ! Cool skin: 161 LOGICAL :: l_use_skin = .FALSE. 162 REAL(wp), DIMENSION(jpi,jpj) :: zsst ! to back up the initial bulk SST 163 164 ! 124 165 REAL(wp), DIMENSION(jpi,jpj) :: func_m, func_h 125 166 REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2 126 167 !!---------------------------------------------------------------------------------- 127 168 ! 169 ! Cool skin ? 170 IF( PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) ) THEN 171 l_use_skin = .TRUE. 172 END IF 173 IF (lwp) PRINT *, ' *** LOLO(sbcblk_algo_ecmwf.F90) => l_use_skin =', l_use_skin 174 128 175 ! Identical first gess as in COARE, with IFS parameter values though 129 176 ! … … 131 178 IF( ABS(zu - zt) < 0.01 ) l_zt_equal_zu = .TRUE. ! testing "zu == zt" is risky with double precision 132 179 180 !! Initialization for cool skin: 181 IF( l_use_skin ) THEN 182 zsst = T_s ! save the bulk SST 183 T_s = T_s - 0.25 ! First guess of correction 184 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 185 END IF 133 186 134 187 !! First guess of temperature and humidity at height zu: 135 t_zu = MAX( t_zt , 0.0 ) ! who knows what's given on masked-continental regions...136 q_zu = MAX( q_zt , 1.e-6 ) ! "188 t_zu = MAX( t_zt , 0.0_wp ) ! who knows what's given on masked-continental regions... 189 q_zu = MAX( q_zt , 1.e-6_wp) ! " 137 190 138 191 !! Pot. temp. difference (and we don't want it to be 0!) 139 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.e-6), dt_zu )140 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.e-9), dq_zu )192 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 193 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 141 194 142 195 znu_a = visc_air(t_zt) ! Air viscosity (m^2/s) at zt given from temperature in (K) 143 196 144 ztmp2 = 0.5 * 0.5! initial guess for wind gustiness contribution197 ztmp2 = 0.5_wp*0.5_wp ! initial guess for wind gustiness contribution 145 198 U_blk = SQRT(U_zu*U_zu + ztmp2) 146 199 147 ! z0 = 0.0001 148 ztmp2 = 10000. ! optimization: ztmp2 == 1/z0 200 ztmp2 = 10000._wp ! optimization: ztmp2 == 1/z0 (with z0 first guess == 0.0001) 149 201 ztmp0 = LOG(zu*ztmp2) 150 202 ztmp1 = LOG(10.*ztmp2) 151 u_star = 0.035*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 152 153 z0 = charn0*u_star*u_star/grav + 0.11*znu_a/u_star 154 z0t = 0.1*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ! WARNING: 1/z0t ! 155 156 Cd = (vkarmn/ztmp0)**2 ! first guess of Cd 157 158 ztmp0 = vkarmn*vkarmn/LOG(zt*z0t)/Cd 203 u_star = 0.035_wp*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 204 205 z0 = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 206 z0 = MIN(ABS(z0), 0.001_wp) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 207 z0t = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) ) 208 z0t = MIN(ABS(z0t), 0.001_wp) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 209 210 ztmp2 = vkarmn/ztmp0 211 Cd = ztmp2*ztmp2 ! first guess of Cd 212 213 ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd 159 214 160 215 ztmp2 = Ri_bulk( zu, t_zu, dt_zu, q_zu, dq_zu, U_blk ) ! Ribu = Bulk Richardson number 161 216 162 217 !! First estimate of zeta_u, depending on the stability, ie sign of Ribu (ztmp2): 163 ztmp1 = 0.5 + SIGN( 0.5 , ztmp2 )218 ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 ) 164 219 func_m = ztmp0*ztmp2 ! temporary array !! 165 !! Ribu < 0 Ribu > 0 Beta = 1.25166 func_h = (1.-ztmp1)*(func_m/(1.+ztmp2/(-zu/(zi0*0.004*Beta0**3)))) & ! temporary array !!! func_h == zeta_u167 & + ztmp1*(func_m*(1. + 27./9.*ztmp2/ztmp0))220 func_h = (1.-ztmp1) * (func_m/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & ! Ribu < 0 ! temporary array !!! func_h == zeta_u 221 & + ztmp1 * (func_m*(1._wp + 27._wp/9._wp*ztmp2/func_m)) ! Ribu > 0 222 !#LB: should make sure that the "func_m" of "27./9.*ztmp2/func_m" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" ! 168 223 169 224 !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L 170 ztmp0 = vkarmn/(LOG(zu*z0t) - psi_h_ecmwf(func_h))225 ztmp0 = vkarmn/(LOG(zu/z0t) - psi_h_ecmwf(func_h)) 171 226 172 227 u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_ecmwf(func_h)) … … 176 231 ! What's need to be done if zt /= zu: 177 232 IF( .NOT. l_zt_equal_zu ) THEN 178 !179 233 !! First update of values at zu (or zt for wind) 180 234 ztmp0 = psi_h_ecmwf(func_h) - psi_h_ecmwf(zt*func_h/zu) ! zt*func_h/zu == zeta_t 181 ztmp1 = log(zt/zu) + ztmp0235 ztmp1 = LOG(zt/zu) + ztmp0 182 236 t_zu = t_zt - t_star/vkarmn*ztmp1 183 237 q_zu = q_zt - q_star/vkarmn*ztmp1 184 q_zu = (0.5 + sign(0.5,q_zu))*q_zu !Makes it impossible to have negative humidity : 185 186 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu ) 187 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu ) 238 q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity : 188 239 ! 189 ENDIF 240 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 241 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 242 END IF 190 243 191 244 … … 195 248 !! First guess of inverse of Monin-Obukov length (1/L) : 196 249 ztmp0 = (1. + rctv0*q_zu) ! the factor to apply to temp. to get virt. temp... 197 Linv = grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / ( u_star*u_star * t_zu*ztmp0 ) 250 Linv = grav*vkarmn*(t_star*ztmp0 + rctv0*t_zu*q_star) / MAX( u_star*u_star * t_zu*ztmp0 , 1.E-9 ) ! #LOLO 251 Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 198 252 199 253 !! Functions such as u* = U_blk*vkarmn/func_m 200 ztmp1 = zu + z0 201 ztmp0 = ztmp1*Linv 202 func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0*Linv) 203 func_h = LOG(ztmp1*z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(1./z0t*Linv) 204 254 ztmp0 = zu*Linv 255 func_m = LOG(zu) - LOG(z0) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf( z0*Linv) 256 func_h = LOG(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 205 257 206 258 !! ITERATION BLOCK 207 !! ***************208 209 259 DO j_itt = 1, nb_itt 210 260 … … 213 263 214 264 !! New estimate of the inverse of the Monin-Obukhon length (Linv == zeta/zu) : 215 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3, p.33, IFS doc - Cy31r1 265 Linv = ztmp0*func_m*func_m/func_h / zu ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 266 !! Note: it is slightly different that the L we would get with the usual 267 !! expression, as in coare algorithm or in 'mod_phymbl.f90' (One_on_L_MO()) 268 Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 216 269 217 270 !! Update func_m with new Linv: 218 ztmp1 = zu + z0 219 func_m = LOG(ztmp1) -LOG(z0) - psi_m_ecmwf(ztmp1*Linv) + psi_m_ecmwf(z0*Linv) 271 func_m = LOG(zu) -LOG(z0) - psi_m_ecmwf(zu*Linv) + psi_m_ecmwf(z0*Linv) ! LB: should be "zu+z0" rather than "zu" alone, but z0 is tiny wrt zu! 220 272 221 273 !! Need to update roughness lengthes: … … 223 275 ztmp2 = u_star*u_star 224 276 ztmp1 = znu_a/u_star 225 z0 = alpha_M*ztmp1 + charn0*ztmp2/grav226 z0t = alpha_H*ztmp1! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1227 z0q = alpha_Q*ztmp1228 277 z0 = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 278 z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 279 z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp) 280 229 281 !! Update wind at 10m taking into acount convection-related wind gustiness: 282 !! => Chap. 3.2, IFS doc - Cy40r1, Eq.3.17 and Eq.3.18 + Eq.3.8 230 283 ! Only true when unstable (L<0) => when ztmp0 < 0 => - !!! 231 ztmp2 = ztmp2 * ( MAX(-zi0*Linv/vkarmn,0.))**(2./3.) ! => w*^2 (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1)284 ztmp2 = ztmp2 * ( MAX(-zi0*Linv/vkarmn , 0._wp))**(2._wp/3._wp) ! => w*^2 (combining Eq. 3.8 and 3.18, hap.3, IFS doc - Cy31r1) 232 285 !! => equivalent using Beta=1 (gustiness parameter, 1.25 for COARE, also zi0=600 in COARE..) 233 U_blk = MAX( sqrt(U_zu*U_zu + ztmp2), 0.2) ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1286 U_blk = MAX( SQRT(U_zu*U_zu + ztmp2) , 0.2_wp ) ! eq.3.17, Chap.3, p.32, IFS doc - Cy31r1 234 287 ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0. 235 288 … … 238 291 !! as well the air-sea differences: 239 292 IF( .NOT. l_zt_equal_zu ) THEN 240 241 293 !! Arrays func_m and func_h are free for a while so using them as temporary arrays... 242 func_h = psi_h_ecmwf( (zu+z0)*Linv) ! temporary array !!!243 func_m = psi_h_ecmwf( (zt+z0)*Linv) ! temporary array !!!294 func_h = psi_h_ecmwf(zu*Linv) ! temporary array !!! 295 func_m = psi_h_ecmwf(zt*Linv) ! temporary array !!! 244 296 245 297 ztmp2 = psi_h_ecmwf(z0t*Linv) 246 298 ztmp0 = func_h - ztmp2 247 ztmp1 = vkarmn/(LOG(zu +z0) - LOG(z0t) - ztmp0)299 ztmp1 = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) 248 300 t_star = dt_zu*ztmp1 249 301 ztmp2 = ztmp0 - func_m + ztmp2 … … 253 305 ztmp2 = psi_h_ecmwf(z0q*Linv) 254 306 ztmp0 = func_h - ztmp2 255 ztmp1 = vkarmn/(LOG(zu +z0) - LOG(z0q) - ztmp0)307 ztmp1 = vkarmn/(LOG(zu) - LOG(z0q) - ztmp0) 256 308 q_star = dq_zu*ztmp1 257 309 ztmp2 = ztmp0 - func_m + ztmp2 258 ztmp1 = log(zt/zu) + ztmp2310 ztmp1 = LOG(zt/zu) + ztmp2 259 311 q_zu = q_zt - q_star/vkarmn*ztmp1 260 261 dt_zu = t_zu - sst ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6), dt_zu )262 dq_zu = q_zu - ssq ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9), dq_zu )263 264 312 END IF 265 313 266 314 !! Updating because of updated z0 and z0t and new Linv... 267 ztmp1 = zu + z0 268 ztmp0 = ztmp1*Linv 269 func_m = log(ztmp1) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 270 func_h = log(ztmp1) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 271 272 END DO 315 ztmp0 = zu*Linv 316 func_m = log(zu) - LOG(z0 ) - psi_m_ecmwf(ztmp0) + psi_m_ecmwf(z0 *Linv) 317 func_h = log(zu) - LOG(z0t) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0t*Linv) 318 319 !! SKIN related part 320 !! ----------------- 321 IF( l_use_skin ) THEN 322 !! compute transfer coefficients at zu : lolo: verifier... 323 Cd = vkarmn*vkarmn/(func_m*func_m) 324 Ch = vkarmn*vkarmn/(func_m*func_h) 325 ztmp1 = LOG(zu) - LOG(z0q) - psi_h_ecmwf(ztmp0) + psi_h_ecmwf(z0q*Linv) ! func_q 326 Ce = vkarmn*vkarmn/(func_m*ztmp1) 327 ! Non-Solar heat flux to the ocean: 328 ztmp1 = U_blk*MAX(rho_air(t_zu, q_zu, slp), 1._wp) ! rho*U10 329 ztmp2 = T_s*T_s 330 ztmp1 = ztmp1 * ( Ce*rLevap*(q_zu - q_s) + Ch*rCp_dry*(t_zu - T_s) ) & ! Total turb. heat flux 331 & + (rad_lw - emiss_w*stefan*ztmp2*ztmp2) ! Net longwave flux 332 !! Updating the values of the skin temperature T_s and q_s : 333 CALL CSWL_ECMWF( Qsw, ztmp1, u_star, zsst, T_s ) 334 q_s = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! 200 -> just to avoid numerics problem on masked regions if silly values are given 335 END IF 336 337 IF( (l_use_skin).OR.(.NOT. l_zt_equal_zu) ) THEN 338 dt_zu = t_zu - T_s ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 339 dq_zu = q_zu - q_s ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 340 END IF 341 342 END DO !DO j_itt = 1, nb_itt 273 343 274 344 Cd = vkarmn*vkarmn/(func_m*func_m) 275 345 Ch = vkarmn*vkarmn/(func_m*func_h) 276 ztmp1 = log((zu + z0)/z0q) - psi_h_ecmwf((zu + z0)*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q 277 Ce = vkarmn*vkarmn/(func_m*ztmp1) 278 279 ztmp1 = zu + z0 280 Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 )) 281 Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t)) 282 Cen = vkarmn*vkarmn / (log(ztmp1/z0q)*log(ztmp1/z0q)) 346 ztmp2 = log(zu/z0q) - psi_h_ecmwf(zu*Linv) + psi_h_ecmwf(z0q*Linv) ! func_q 347 Ce = vkarmn*vkarmn/(func_m*ztmp2) 348 349 Cdn = vkarmn*vkarmn / (log(zu/z0 )*log(zu/z0 )) 350 Chn = vkarmn*vkarmn / (log(zu/z0t)*log(zu/z0t)) 351 Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 283 352 284 353 END SUBROUTINE TURB_ECMWF … … 294 363 !! and L is M-O length 295 364 !! 296 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)365 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 297 366 !!---------------------------------------------------------------------------------- 298 367 REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ecmwf … … 306 375 DO ji = 1, jpi 307 376 ! 308 zzeta = MIN( pzeta(ji,jj) , 5. ) !! Very stable conditions (L positif and big!):377 zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 309 378 ! 310 379 ! Unstable (Paulson 1970): 311 380 ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 312 zx = SQRT(ABS(1. - 16.*zzeta))313 ztmp = 1. + SQRT(zx)381 zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 382 ztmp = 1._wp + SQRT(zx) 314 383 ztmp = ztmp*ztmp 315 psi_unst = LOG( 0.125 *ztmp*(1.+ zx) ) &316 & -2. *ATAN( SQRT(zx) ) + 0.5*rpi384 psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) ) & 385 & -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 317 386 ! 318 387 ! Unstable: 319 388 ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 320 psi_stab = -2. /3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) &321 & - zzeta - 2. /3.*5./0.35389 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 390 & - zzeta - 2._wp/3._wp*5._wp/0.35_wp 322 391 ! 323 392 ! Combining: 324 stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1325 ! 326 psi_m_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable327 & + stab * psi_stab ! (zzeta > 0) Stable393 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 394 ! 395 psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 396 & + stab * psi_stab ! (zzeta > 0) Stable 328 397 ! 329 398 END DO … … 332 401 END FUNCTION psi_m_ecmwf 333 402 334 403 335 404 FUNCTION psi_h_ecmwf( pzeta ) 336 405 !!---------------------------------------------------------------------------------- … … 342 411 !! and L is M-O length 343 412 !! 344 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)413 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 345 414 !!---------------------------------------------------------------------------------- 346 415 REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ecmwf … … 354 423 DO ji = 1, jpi 355 424 ! 356 zzeta = MIN(pzeta(ji,jj) , 5. ) ! Very stable conditions (L positif and big!):357 ! 358 zx = ABS(1. - 16.*zzeta)**.25 ! this is actually (1/phi_m)**2 !!!425 zzeta = MIN(pzeta(ji,jj) , 5._wp) ! Very stable conditions (L positif and big!): 426 ! 427 zx = ABS(1._wp - 16._wp*zzeta)**.25 ! this is actually (1/phi_m)**2 !!! 359 428 ! ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 360 429 ! Unstable (Paulson 1970) : 361 psi_unst = 2. *LOG(0.5*(1.+ zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1430 psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx)) ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 362 431 ! 363 432 ! Stable: 364 psi_stab = -2. /3.*(zzeta - 5./0.35)*EXP(-0.35*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1365 & - ABS(1. + 2./3.*zzeta)**1.5 - 2./3.*5./0.35 + 1.433 psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 434 & - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 366 435 ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 367 436 ! 368 stab = 0.5 + SIGN(0.5, zzeta) ! zzeta > 0 => stab = 1369 ! 370 ! 371 psi_h_ecmwf(ji,jj) = (1. - stab) * psi_unst & ! (zzeta < 0) Unstable372 & + stab * psi_stab ! (zzeta > 0) Stable437 stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 438 ! 439 ! 440 psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 441 & + stab * psi_stab ! (zzeta > 0) Stable 373 442 ! 374 443 END DO … … 382 451 !! Bulk Richardson number (Eq. 3.25 IFS doc) 383 452 !! 384 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)453 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 385 454 !!---------------------------------------------------------------------------------- 386 455 REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk ! … … 394 463 !!---------------------------------------------------------------------------------- 395 464 ! 396 Ri_bulk = grav*pz/(pub*pub) 397 & * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/( Cp_dry+Cp_vap*pqz)))&465 Ri_bulk = grav*pz/(pub*pub) & 466 & * ( pdt/(ptz - 0.5_wp*(pdt + grav*pz/(rCp_dry + rCp_vap*pqz))) & 398 467 & + rctv0*pdq ) 399 468 ! … … 401 470 402 471 403 FUNCTION visc_air(ptak)404 !!----------------------------------------------------------------------------------405 !! Air kinetic viscosity (m^2/s) given from temperature in degrees...406 !!407 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)408 !!----------------------------------------------------------------------------------409 REAL(wp), DIMENSION(jpi,jpj) :: visc_air !410 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K)411 !412 INTEGER :: ji, jj ! dummy loop indices413 REAL(wp) :: ztc, ztc2 ! local scalar414 !!----------------------------------------------------------------------------------415 !416 DO jj = 1, jpj417 DO ji = 1, jpi418 ztc = ptak(ji,jj) - rt0 ! air temp, in deg. C419 ztc2 = ztc*ztc420 visc_air(ji,jj) = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)421 END DO422 END DO423 !424 END FUNCTION visc_air425 472 426 473 !!====================================================================== -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90
r10190 r11111 11 11 !! 12 12 !! Routine turb_ncar maintained and developed in AeroBulk 13 !! (http ://aerobulk.sourceforge.net/)13 !! (https://github.com/brodeau/aerobulk/) 14 14 !! 15 15 !! L. Brodeau, 2015 … … 43 43 44 44 PUBLIC :: TURB_NCAR ! called by sbcblk.F90 45 46 ! ! NCAR own values for given constants:47 REAL(wp), PARAMETER :: rctv0 = 0.608 ! constant to obtain virtual temperature...48 45 49 46 !!---------------------------------------------------------------------- … … 76 73 !! => 'vkarmn' and 'grav' 77 74 !! 78 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)75 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 79 76 !! 80 77 !! INPUT : … … 238 235 !! Origin: Large & Yeager 2008 eq.(11a) and eq.(11b) 239 236 !! 240 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)237 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 241 238 !!---------------------------------------------------------------------------------- 242 239 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pw10 ! scalar wind speed at 10m (m/s) … … 277 274 !! and L is M-O length 278 275 !! 279 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)276 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 280 277 !!---------------------------------------------------------------------------------- 281 278 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta … … 311 308 !! and L is M-O length 312 309 !! 313 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https:// sourceforge.net/p/aerobulk)310 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 314 311 !!---------------------------------------------------------------------------------- 315 312 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
Note: See TracChangeset
for help on using the changeset viewer.
