- Timestamp:
- 2019-07-18T11:35:28+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC
- Files:
-
- 2 added
- 2 deleted
- 2 edited
-
sbcblk.F90 (modified) (17 diffs)
-
sbcblk_algo_coare.F90 (deleted)
-
sbcblk_algo_coare3p0.F90 (added)
-
sbcblk_algo_coare3p5.F90 (deleted)
-
sbcblk_algo_coare3p6.F90 (added)
-
sbcblk_phy.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90
r11266 r11284 46 46 #endif 47 47 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 48 USE sbcblk_algo_coare ! => turb_coare: COAREv3.0 (Fairall et al. 2003)49 USE sbcblk_algo_coare3p 5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)48 USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003) 49 USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013) 50 50 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31) 51 51 ! … … 87 87 LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008) 88 88 LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003) 89 LOGICAL :: ln_COARE_3p 5 ! "COARE 3.5" algorithm (Edson et al. 2013)89 LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013) 90 90 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 31) 91 91 ! … … 109 109 !LB: 110 110 LOGICAL :: ln_skin ! use the cool-skin / warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 111 LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature if .true. (specific humidity espected if .false.) !LB 111 LOGICAL :: ln_humi_sph ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB 112 LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 113 LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB 114 ! 115 INTEGER :: nhumi ! choice of the bulk algorithm 116 ! ! associated indices: 117 INTEGER, PARAMETER :: np_humi_sph = 1 118 INTEGER, PARAMETER :: np_humi_dpt = 2 119 INTEGER, PARAMETER :: np_humi_rlh = 3 112 120 !LB. 113 121 … … 116 124 INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008) 117 125 INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003) 118 INTEGER, PARAMETER :: np_COARE_3p 5 = 3 ! "COARE 3.5" algorithm (Edson et al. 2013)126 INTEGER, PARAMETER :: np_COARE_3p6 = 3 ! "COARE 3.6" algorithm (Edson et al. 2013) 119 127 INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 31) 120 128 … … 172 180 NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields 173 181 & sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, & 174 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p 5, ln_ECMWF, & ! bulk algorithm182 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm 175 183 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 176 184 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, & 177 & ln_skin, ln_humi_ dpt ! cool-skin / warm-layer param.!LB185 & ln_skin, ln_humi_sph, ln_humi_dpt, ln_humi_rlh ! cool-skin / warm-layer !LB 178 186 !!--------------------------------------------------------------------- 179 187 ! … … 201 209 nblk = np_COARE_3p0 ; ioptio = ioptio + 1 202 210 ENDIF 203 IF( ln_COARE_3p 5) THEN204 nblk = np_COARE_3p 5; ioptio = ioptio + 1211 IF( ln_COARE_3p6 ) THEN 212 nblk = np_COARE_3p6 ; ioptio = ioptio + 1 205 213 ENDIF 206 214 IF( ln_ECMWF ) THEN 207 215 nblk = np_ECMWF ; ioptio = ioptio + 1 208 216 ENDIF 209 ! 217 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 218 219 220 210 221 !LB: 211 222 ! !** initialization of the cool-skin / warm-layer parametrization … … 215 226 IF( sbc_blk_cswl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 216 227 END IF 228 ! 229 ioptio = 0 230 IF( ln_humi_sph ) THEN 231 nhumi = np_humi_sph ; ioptio = ioptio + 1 232 ENDIF 233 IF( ln_humi_dpt ) THEN 234 nhumi = np_humi_dpt ; ioptio = ioptio + 1 235 ENDIF 236 IF( ln_humi_rlh ) THEN 237 nhumi = np_humi_rlh ; ioptio = ioptio + 1 238 ENDIF 239 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' ) 217 240 !LB. 218 241 219 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )220 242 ! 221 243 IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr … … 278 300 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR 279 301 WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0 280 WRITE(numout,*) ' "COARE 3. 5" algorithm (Edson et al. 2013) ln_COARE_3p5 = ', ln_COARE_3p5302 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 281 303 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 31) ln_ECMWF = ', ln_ECMWF 282 304 WRITE(numout,*) ' add High freq.contribution to the stress module ln_taudif = ', ln_taudif … … 294 316 CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)' 295 317 CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)' 296 CASE( np_COARE_3p 5 ) ; WRITE(numout,*) ' ==>>> "COARE 3.5" algorithm (Edson et al. 2013)'318 CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 297 319 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 31)' 298 320 END SELECT … … 300 322 WRITE(numout,*) 301 323 WRITE(numout,*) ' use cool-skin/warm-layer parameterization (SSST) ln_skin = ', ln_skin !LB 324 ! 325 !LB: 302 326 WRITE(numout,*) 303 WRITE(numout,*) ' air humidity read in files is dew-point temperature? ln_humi_dpt = ', ln_humi_dpt !LB 327 SELECT CASE( nhumi ) !* Print the choice of air humidity 328 CASE( np_humi_sph ) ; WRITE(numout,*) ' ==>>> air humidity is SPECIFIC HUMIDITY [kg/kg]' 329 CASE( np_humi_dpt ) ; WRITE(numout,*) ' ==>>> air humidity is DEW-POINT TEMPERATURE [K]' 330 CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]' 331 END SELECT 332 !LB. 304 333 ! 305 334 ENDIF … … 360 389 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 361 390 !LB: 362 IF ( ln_humi_dpt ) THEN 363 qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 364 ELSE 365 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 366 END IF 391 SELECT CASE( nhumi ) 392 CASE( np_humi_sph ) 393 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 394 CASE( np_humi_dpt ) 395 qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 396 CASE( np_humi_rlh ) 397 qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !RH is % percent in file 398 END SELECT 367 399 !LB. 368 400 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac … … 471 503 472 504 !LB: 473 ! spec. hum. of air at "rn_zqt" m above thes sea: 474 IF ( ln_humi_dpt ) THEN 475 ! spec. hum. of air is deduced from dew-point temperature and SLP: q_air = q_sat( d_air, slp) 505 ! zqair = specific humidity of air at "rn_zqt" m above thes sea: 506 SELECT CASE( nhumi ) 507 CASE( np_humi_sph ) 508 zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity! 509 CASE( np_humi_dpt ) 476 510 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of d_air and slp !' 477 511 zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 478 ELSE 479 zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity! 480 END IF 512 CASE( np_humi_rlh ) 513 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of RH, t_air and slp !' 514 zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 515 END SELECT 481 516 !LB. 482 517 … … 485 520 !! (since reanalysis products provide T at z, not theta !) 486 521 ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt 487 522 488 523 489 524 IF ( ln_skin ) THEN 490 525 491 526 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 492 527 493 528 CASE( np_COARE_3p0 ) 494 CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.0 529 CALL turb_coare3p0 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.0 530 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 531 & Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & 532 & Tsk_b=tsk(:,:) ) 533 534 CASE( np_COARE_3p6 ) 535 CALL turb_coare3p6 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.6 495 536 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 496 537 & Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & … … 502 543 & Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), & 503 544 & Tsk_b=tsk(:,:) ) 504 545 505 546 CASE DEFAULT 506 547 CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin==.true."' ) … … 523 564 ELSE !IF ( ln_skin ) 524 565 525 566 526 567 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 527 568 ! … … 530 571 531 572 CASE( np_COARE_3p0 ) 532 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare " WITHOUT CSWL optional arrays!!!'533 CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.0573 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITHOUT CSWL optional arrays!!!' 574 CALL turb_coare3p0 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.0 534 575 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 535 576 536 CASE( np_COARE_3p 5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.5577 CASE( np_COARE_3p6 ) ; CALL turb_coare3p6( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, & ! COARE v3.6 537 578 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 538 579 … … 802 843 803 844 !LB: 804 ! spec. hum. of air at "rn_zqt" m above thes sea: 805 IF ( ln_humi_dpt ) THEN 806 ! spec. hum. of air is deduced from dew-point temperature and SLP: q_air = q_sat( d_air, slp) 845 SELECT CASE( nhumi ) 846 CASE( np_humi_sph ) 847 zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity! 848 CASE( np_humi_dpt ) 807 849 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of d_air and slp !' 808 850 zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 809 ELSE 810 zqair(:,:) = sf(jp_humi)%fnow(:,:,1) ! what we read in file is already a spec. humidity! 811 END IF 851 CASE( np_humi_rlh ) 852 IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of RH, t_air and slp !' 853 zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 854 END SELECT 812 855 !LB. 813 856 -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90
r11266 r11284 17 17 !! Ri_bulk : bulk Richardson number aka BRN 18 18 !! q_sat : saturation humidity as a function of SLP and temperature 19 19 !! q_air_rh : specific humidity as a function of RH, t_air and SLP 20 20 21 USE dom_oce ! ocean space and time domain 21 22 USE phycst ! physical constants … … 37 38 PUBLIC Ri_bulk 38 39 PUBLIC q_sat 40 PUBLIC q_air_rh 39 41 40 42 !!---------------------------------------------------------------------- … … 115 117 !! *** FUNCTION L_vap *** 116 118 !! 117 !! ** Purpose : Compute the latent heat of vaporization of water fromtemperature119 !! ** Purpose : Compute the latent heat of vaporization of water out of temperature 118 120 !! 119 121 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) … … 273 275 274 276 277 FUNCTION e_sat_sclr( ptak, pslp ) 278 !!---------------------------------------------------------------------------------- 279 !! *** FUNCTION e_sat_sclr *** 280 !! < SCALAR argument version > 281 !! ** Purpose : water vapor at saturation in [Pa] 282 !! Based on accurate estimate by Goff, 1957 283 !! 284 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 285 !! 286 !! Note: what rt0 should be here, is 273.16 (triple point of water) and not 273.15 like here 287 !!---------------------------------------------------------------------------------- 288 REAL(wp), INTENT(in) :: ptak ! air temperature [K] 289 REAL(wp), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa] 290 REAL(wp) :: e_sat_sclr ! water vapor at saturation [kg/kg] 291 ! 292 REAL(wp) :: zta, ztmp ! local scalar 293 !!---------------------------------------------------------------------------------- 294 ! 295 zta = MAX( ptak , 180._wp ) ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions... 296 ztmp = rt0 / zta 297 ! 298 ! Vapour pressure at saturation [Pa] : WMO, (Goff, 1957) 299 e_sat_sclr = 100.*( 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0) & 300 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) ) & 301 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) ) 302 ! 303 END FUNCTION e_sat_sclr 304 305 275 306 FUNCTION q_sat( ptak, pslp ) 276 307 !!---------------------------------------------------------------------------------- … … 288 319 ! 289 320 INTEGER :: ji, jj ! dummy loop indices 290 REAL(wp) :: zta, ze_sat, ztmp ! local scalar 291 !!---------------------------------------------------------------------------------- 292 ! 293 DO jj = 1, jpj 294 DO ji = 1, jpi 295 ! 296 zta = MAX( ptak(ji,jj) , 180._wp ) ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions... 297 ztmp = rt0 / zta 298 ! 299 ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 300 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0) & 301 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) ) & 302 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 ) 303 ! 304 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] 321 REAL(wp) :: ze_sat ! local scalar 322 !!---------------------------------------------------------------------------------- 323 ! 324 DO jj = 1, jpj 325 DO ji = 1, jpi 326 ! 327 ze_sat = e_sat_sclr( ptak(ji,jj), pslp(ji,jj) ) 328 ! 329 q_sat(ji,jj) = reps0 * ze_sat/( pslp(ji,jj) - (1._wp - reps0)*ze_sat ) 305 330 ! 306 331 END DO … … 308 333 ! 309 334 END FUNCTION q_sat 335 336 FUNCTION q_air_rh(prha, ptak, pslp) 337 !!---------------------------------------------------------------------------------- 338 !! Specific humidity of air out of Relative Humidity 339 !! 340 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 341 !!---------------------------------------------------------------------------------- 342 REAL(wp), DIMENSION(jpi,jpj) :: q_air_rh 343 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prha !: relative humidity [fraction, not %!!!] 344 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak !: air temperature [K] 345 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp !: atmospheric pressure [Pa] 346 ! 347 INTEGER :: ji, jj ! dummy loop indices 348 REAL(wp) :: ze ! local scalar 349 !!---------------------------------------------------------------------------------- 350 ! 351 DO jj = 1, jpj 352 DO ji = 1, jpi 353 ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj), pslp(ji,jj)) 354 q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze) 355 END DO 356 END DO 357 ! 358 END FUNCTION q_air_rh 310 359 311 360
Note: See TracChangeset
for help on using the changeset viewer.
