- Timestamp:
- 2021-03-05T17:03:57+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/SBC/sbc_phy.F90
r14110 r14592 23 23 24 24 IMPLICIT NONE 25 !PRIVATE 26 PUBLIC !! Haleluja that was the solution 25 PUBLIC !! Haleluja that was the solution for AGRIF 27 26 28 27 INTEGER , PARAMETER, PUBLIC :: nb_iter0 = 5 ! Default number of itterations in bulk-param algorithms (can be overriden b.m.o `nb_iter` optional argument) 29 28 30 29 !! (mainly removed from sbcblk.F90) 31 REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp !: Specic heat of dry air, constant pressure [J/K/kg] 32 REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp !: Specic heat of water vapor, constant pressure [J/K/kg] 33 REAL(wp), PARAMETER, PUBLIC :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg] 34 REAL(wp), PARAMETER, PUBLIC :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg] 35 REAL(wp), PARAMETER, PUBLIC :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622 36 REAL(wp), PARAMETER, PUBLIC :: rctv0 = R_vap/R_dry - 1._wp !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 37 REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp !: specific heat of air (only used for ice fluxes now...) 38 REAL(wp), PARAMETER, PUBLIC :: albo = 0.066_wp !: ocean albedo assumed to be constant 30 REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp !: Specic heat of dry air, constant pressure [J/K/kg] 31 REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp !: Specic heat of water vapor, constant pressure [J/K/kg] 32 REAL(wp), PARAMETER, PUBLIC :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg] 33 REAL(wp), PARAMETER, PUBLIC :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg] 34 REAL(wp), PARAMETER, PUBLIC :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622 35 REAL(wp), PARAMETER, PUBLIC :: rctv0 = R_vap/R_dry - 1._wp !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 36 REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp !: specific heat of air (only used for ice fluxes now...) 37 REAL(wp), PARAMETER, PUBLIC :: albo = 0.066_wp !: ocean albedo assumed to be constant 38 REAL(wp), PARAMETER, PUBLIC :: R_gas = 8.314510_wp !: Universal molar gas constant [J/mol/K] 39 REAL(wp), PARAMETER, PUBLIC :: rmm_dryair = 28.9647e-3_wp !: dry air molar mass / molecular weight [kg/mol] 40 REAL(wp), PARAMETER, PUBLIC :: rmm_water = 18.0153e-3_wp !: water molar mass / molecular weight [kg/mol] 41 REAL(wp), PARAMETER, PUBLIC :: rmm_ratio = rmm_water / rmm_dryair 42 REAL(wp), PARAMETER, PUBLIC :: rgamma_dry = R_gas / ( rmm_dryair * rCp_dry ) !: Poisson constant for dry air 43 REAL(wp), PARAMETER, PUBLIC :: rpref = 1.e5_wp !: reference air pressure for exner function [Pa] 39 44 ! 40 45 REAL(wp), PARAMETER, PUBLIC :: rho0_a = 1.2_wp !: Approx. of density of air [kg/m^3] … … 83 88 END INTERFACE virt_temp 84 89 90 INTERFACE pres_temp 91 MODULE PROCEDURE pres_temp_vctr, pres_temp_sclr 92 END INTERFACE pres_temp 93 94 INTERFACE theta_exner 95 MODULE PROCEDURE theta_exner_vctr, theta_exner_sclr 96 END INTERFACE theta_exner 97 85 98 INTERFACE visc_air 86 99 MODULE PROCEDURE visc_air_vctr, visc_air_sclr … … 98 111 MODULE PROCEDURE e_sat_ice_vctr, e_sat_ice_sclr 99 112 END INTERFACE e_sat_ice 113 100 114 INTERFACE de_sat_dt_ice 101 115 MODULE PROCEDURE de_sat_dt_ice_vctr, de_sat_dt_ice_sclr … … 148 162 149 163 PUBLIC virt_temp 164 PUBLIC pres_temp 165 PUBLIC theta_exner 150 166 PUBLIC rho_air 151 167 PUBLIC visc_air … … 158 174 PUBLIC q_air_rh 159 175 PUBLIC dq_sat_dt_ice 160 ! :176 ! 161 177 PUBLIC update_qnsol_tau 162 178 PUBLIC alpha_sw … … 205 221 ! 206 222 END FUNCTION virt_temp_sclr 207 !! 223 208 224 FUNCTION virt_temp_vctr( pta, pqa ) 225 209 226 REAL(wp), DIMENSION(jpi,jpj) :: virt_temp_vctr !: virtual temperature [K] 210 227 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K] 211 228 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air [kg/kg] 229 212 230 virt_temp_vctr(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:)) 231 213 232 END FUNCTION virt_temp_vctr 214 !=============================================================================================== 233 234 235 FUNCTION pres_temp_sclr( pqspe, pslp, pz, ptpot, pta, l_ice ) 236 237 !!------------------------------------------------------------------------------- 238 !! *** FUNCTION pres_temp *** 239 !! 240 !! ** Purpose : compute air pressure using barometric equation 241 !! from either potential or absolute air temperature 242 !! ** Author: G. Samson, Feb 2021 243 !!------------------------------------------------------------------------------- 244 245 REAL(wp) :: pres_temp_sclr ! air pressure [Pa] 246 REAL(wp), INTENT(in ) :: pqspe ! air specific humidity [kg/kg] 247 REAL(wp), INTENT(in ) :: pslp ! sea-level pressure [Pa] 248 REAL(wp), INTENT(in ) :: pz ! height above surface [m] 249 REAL(wp), INTENT(in ) , OPTIONAL :: ptpot ! air potential temperature [K] 250 REAL(wp), INTENT(inout), OPTIONAL :: pta ! air absolute temperature [K] 251 REAL(wp) :: ztpot, zta, zpa, zxm, zmask, zqsat 252 INTEGER :: it, niter = 3 ! iteration indice and number 253 LOGICAL , INTENT(in) , OPTIONAL :: l_ice ! sea-ice presence 254 LOGICAL :: lice ! sea-ice presence 255 256 IF( PRESENT(ptpot) ) THEN 257 zmask = 1._wp 258 ztpot = ptpot 259 ELSE 260 zmask = 0._wp 261 zta = pta 262 ENDIF 263 264 lice = .FALSE. 265 IF( PRESENT(l_ice) ) lice = l_ice 266 267 zpa = pslp ! air pressure first guess [Pa] 268 DO it = 1, niter 269 zta = ztpot * ( zpa / rpref )**rgamma_dry * zmask + (1._wp - zmask) * zta 270 zqsat = q_sat( zta, zpa, l_ice=lice ) ! saturation specific humidity [kg/kg] 271 zxm = (1._wp - pqspe/zqsat) * rmm_dryair + pqspe/zqsat * rmm_water ! moist air molar mass [kg/mol] 272 zpa = pslp * EXP( -grav * zxm * pz / ( R_gas * zta ) ) 273 END DO 274 275 pres_temp_sclr = zpa 276 IF(( PRESENT(pta) ).AND.( PRESENT(ptpot) )) pta = zta 277 278 END FUNCTION pres_temp_sclr 279 280 281 FUNCTION pres_temp_vctr( pqspe, pslp, pz, ptpot, pta, l_ice ) 282 283 !!------------------------------------------------------------------------------- 284 !! *** FUNCTION pres_temp *** 285 !! 286 !! ** Purpose : compute air pressure using barometric equation 287 !! from either potential or absolute air temperature 288 !! ** Author: G. Samson, Feb 2021 289 !!------------------------------------------------------------------------------- 290 291 REAL(wp), DIMENSION(jpi,jpj) :: pres_temp_vctr ! air pressure [Pa] 292 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqspe ! air specific humidity [kg/kg] 293 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pslp ! sea-level pressure [Pa] 294 REAL(wp), INTENT(in ) :: pz ! height above surface [m] 295 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) , OPTIONAL :: ptpot ! air potential temperature [K] 296 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout), OPTIONAL :: pta ! air absolute temperature [K] 297 INTEGER :: ji, jj ! loop indices 298 LOGICAL , INTENT(in) , OPTIONAL :: l_ice ! sea-ice presence 299 LOGICAL :: lice ! sea-ice presence 300 301 lice = .FALSE. 302 IF( PRESENT(l_ice) ) lice = l_ice 303 304 IF( PRESENT(ptpot) ) THEN 305 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 306 pres_temp_vctr(ji,jj) = pres_temp_sclr( pqspe(ji,jj), pslp(ji,jj), pz, ptpot=ptpot(ji,jj), pta=pta(ji,jj), l_ice=lice ) 307 END_2D 308 ELSE 309 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 310 pres_temp_vctr(ji,jj) = pres_temp_sclr( pqspe(ji,jj), pslp(ji,jj), pz, pta=pta(ji,jj), l_ice=lice ) 311 END_2D 312 ENDIF 313 314 END FUNCTION pres_temp_vctr 315 316 317 FUNCTION theta_exner_sclr( pta, ppa ) 318 319 !!------------------------------------------------------------------------------- 320 !! *** FUNCTION theta_exner *** 321 !! 322 !! ** Purpose : compute air/surface potential temperature from absolute temperature 323 !! and pressure using Exner function 324 !! ** Author: G. Samson, Feb 2021 325 !!------------------------------------------------------------------------------- 326 327 REAL(wp) :: theta_exner_sclr ! air/surface potential temperature [K] 328 REAL(wp), INTENT(in) :: pta ! air/surface absolute temperature [K] 329 REAL(wp), INTENT(in) :: ppa ! air/surface pressure [Pa] 330 331 theta_exner_sclr = pta * ( rpref / ppa ) ** rgamma_dry 332 333 END FUNCTION theta_exner_sclr 334 335 FUNCTION theta_exner_vctr( pta, ppa ) 336 337 !!------------------------------------------------------------------------------- 338 !! *** FUNCTION theta_exner *** 339 !! 340 !! ** Purpose : compute air/surface potential temperature from absolute temperature 341 !! and pressure using Exner function 342 !! ** Author: G. Samson, Feb 2021 343 !!------------------------------------------------------------------------------- 344 345 REAL(wp), DIMENSION(jpi,jpj) :: theta_exner_vctr ! air/surface potential temperature [K] 346 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta ! air/surface absolute temperature [K] 347 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! air/surface pressure [Pa] 348 INTEGER :: ji, jj ! loop indices 349 350 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 351 theta_exner_vctr(ji,jj) = theta_exner_sclr( pta(ji,jj), ppa(ji,jj) ) 352 END_2D 353 354 END FUNCTION theta_exner_vctr 215 355 216 356 … … 228 368 REAL(wp), DIMENSION(jpi,jpj) :: rho_air_vctr ! density of moist air [kg/m^3] 229 369 !!------------------------------------------------------------------------------- 370 230 371 rho_air_vctr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp ) 372 231 373 END FUNCTION rho_air_vctr 232 374 … … 245 387 !!------------------------------------------------------------------------------- 246 388 rho_air_sclr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp ) 389 247 390 END FUNCTION rho_air_sclr 248 249 250 391 251 392 … … 269 410 270 411 FUNCTION visc_air_vctr(ptak) 412 271 413 REAL(wp), DIMENSION(jpi,jpj) :: visc_air_vctr ! kinetic viscosity (m^2/s) 272 414 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K) 273 415 INTEGER :: ji, jj ! dummy loop indices 274 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 275 visc_air_vctr(ji,jj) = visc_air_sclr( ptak(ji,jj) ) 276 END_2D 416 417 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 418 visc_air_vctr(ji,jj) = visc_air_sclr( ptak(ji,jj) ) 419 END_2D 420 277 421 END FUNCTION visc_air_vctr 278 422 … … 319 463 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 320 464 !!------------------------------------------------------------------------------- 321 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]465 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] 322 466 REAL(wp), DIMENSION(jpi,jpj) :: cp_air_vctr ! specific heat of moist air [J/K/kg] 323 467 !!------------------------------------------------------------------------------- 468 324 469 cp_air_vctr = rCp_dry + rCp_vap * pqa 470 325 471 END FUNCTION cp_air_vctr 326 472 … … 336 482 REAL(wp) :: cp_air_sclr ! specific heat of moist air [J/K/kg] 337 483 !!------------------------------------------------------------------------------- 484 338 485 cp_air_sclr = rCp_dry + rCp_vap * pqa 486 339 487 END FUNCTION cp_air_sclr 340 488 341 489 342 343 !===============================================================================================344 490 FUNCTION gamma_moist_sclr( ptak, pqa ) 345 491 !!---------------------------------------------------------------------------------- … … 366 512 !! 367 513 END FUNCTION gamma_moist_sclr 368 !! 514 369 515 FUNCTION gamma_moist_vctr( ptak, pqa ) 516 370 517 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist_vctr 371 518 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak 372 519 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa 373 520 INTEGER :: ji, jj 374 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 375 gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) ) 376 END_2D 521 522 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 523 gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) ) 524 END_2D 525 377 526 END FUNCTION gamma_moist_vctr 378 !===============================================================================================379 527 380 528 … … 399 547 ! 400 548 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 401 ! 402 zqa = (1._wp + rctv0*pqa(ji,jj)) 403 ! 404 ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with: 405 ! a/ -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one! 406 ! or 407 ! b/ -u* [ theta* + 0.61 theta q* ] 408 ! 409 One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) & 410 & / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) 411 ! 549 zqa = (1._wp + rctv0*pqa(ji,jj)) 550 ! 551 ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with: 552 ! a/ -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one! 553 ! or 554 ! b/ -u* [ theta* + 0.61 theta q* ] 555 ! 556 One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) & 557 & / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp ) 412 558 END_2D 413 559 ! … … 417 563 418 564 419 !===============================================================================================420 565 FUNCTION Ri_bulk_sclr( pz, psst, ptha, pssq, pqa, pub, pta_layer, pqa_layer ) 421 566 !!---------------------------------------------------------------------------------- … … 428 573 REAL(wp) :: Ri_bulk_sclr 429 574 REAL(wp), INTENT(in) :: pz ! height above the sea (aka "delta z") [m] 430 REAL(wp), INTENT(in) :: psst ! SST[K]575 REAL(wp), INTENT(in) :: psst ! potential SST [K] 431 576 REAL(wp), INTENT(in) :: ptha ! pot. air temp. at height "pz" [K] 432 577 REAL(wp), INTENT(in) :: pssq ! 0.98*q_sat(SST) [kg/kg] … … 438 583 LOGICAL :: l_ptqa_l_prvd = .FALSE. 439 584 REAL(wp) :: zqa, zta, zgamma, zdthv, ztv, zsstv ! local scalars 585 REAL(wp) :: ztptv 440 586 !!------------------------------------------------------------------- 441 587 IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE. 442 588 ! 443 zsstv = virt_temp_sclr( psst, pssq ) ! virtual SST (absolute==potential because z=0!) 444 ! 445 zdthv = virt_temp_sclr( ptha, pqa ) - zsstv ! air-sea delta of "virtual potential temperature" 446 ! 447 !! ztv: estimate of the ABSOLUTE virtual temp. within the layer 448 IF( l_ptqa_l_prvd ) THEN 449 ztv = virt_temp_sclr( pta_layer, pqa_layer ) 450 ELSE 451 zqa = 0.5_wp*( pqa + pssq ) ! ~ mean q within the layer... 452 zta = 0.5_wp*( psst + ptha - gamma_moist(ptha, zqa)*pz ) ! ~ mean absolute temperature of air within the layer 453 zta = 0.5_wp*( psst + ptha - gamma_moist( zta, zqa)*pz ) ! ~ mean absolute temperature of air within the layer 454 zgamma = gamma_moist(zta, zqa) ! Adiabatic lapse-rate for moist air within the layer 455 ztv = 0.5_wp*( zsstv + virt_temp_sclr( ptha-zgamma*pz, pqa ) ) 456 END IF 457 ! 458 Ri_bulk_sclr = grav*zdthv*pz / ( ztv*pub*pub ) ! the usual definition of Ri_bulk_sclr 589 zsstv = virt_temp_sclr( psst, pssq ) ! virtual potential SST 590 ztptv = virt_temp_sclr( ptha, pqa ) ! virtual potential air temperature 591 zdthv = ztptv - zsstv ! air-sea delta of "virtual potential temperature" 592 ! 593 Ri_bulk_sclr = grav * zdthv * pz / ( ztptv * pub * pub ) ! the usual definition of Ri_bulk_sclr 459 594 ! 460 595 END FUNCTION Ri_bulk_sclr 461 !! 596 462 597 FUNCTION Ri_bulk_vctr( pz, psst, ptha, pssq, pqa, pub, pta_layer, pqa_layer ) 598 463 599 REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk_vctr 464 600 REAL(wp) , INTENT(in) :: pz ! height above the sea (aka "delta z") [m] … … 473 609 LOGICAL :: l_ptqa_l_prvd = .FALSE. 474 610 INTEGER :: ji, jj 611 475 612 IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE. 476 613 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 477 IF( l_ptqa_l_prvd ) THEN 478 Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), & 479 & pta_layer=pta_layer(ji,jj ), pqa_layer=pqa_layer(ji,jj ) ) 480 ELSE 481 Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) ) 482 END IF 483 END_2D 614 IF( l_ptqa_l_prvd ) THEN !!GS: "IF" inside loop needs to be removed 615 Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), & 616 & pta_layer=pta_layer(ji,jj ), pqa_layer=pqa_layer(ji,jj ) ) 617 ELSE 618 Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) ) 619 END IF 620 END_2D 621 484 622 END FUNCTION Ri_bulk_vctr 485 !=============================================================================================== 486 487 !=============================================================================================== 623 624 488 625 FUNCTION e_sat_sclr( ptak ) 489 626 !!---------------------------------------------------------------------------------- … … 510 647 ! 511 648 END FUNCTION e_sat_sclr 512 !! 649 513 650 FUNCTION e_sat_vctr(ptak) 514 651 REAL(wp), DIMENSION(jpi,jpj) :: e_sat_vctr !: vapour pressure at saturation [Pa] … … 519 656 END_2D 520 657 END FUNCTION e_sat_vctr 521 !=============================================================================================== 522 523 !=============================================================================================== 658 659 524 660 FUNCTION e_sat_ice_sclr(ptak) 525 661 !!--------------------------------------------------------------------------------- … … 537 673 !! 538 674 e_sat_ice_sclr = 100._wp * 10._wp**zle 675 539 676 END FUNCTION e_sat_ice_sclr 540 !! 677 541 678 FUNCTION e_sat_ice_vctr(ptak) 542 679 !! Same as "e_sat" but over ice rather than water! … … 546 683 !!---------------------------------------------------------------------------------- 547 684 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 548 e_sat_ice_vctr(ji,jj) = e_sat_ice_sclr( ptak(ji,jj) ) 549 END_2D 685 e_sat_ice_vctr(ji,jj) = e_sat_ice_sclr( ptak(ji,jj) ) 686 END_2D 687 550 688 END FUNCTION e_sat_ice_vctr 551 !! 689 690 552 691 FUNCTION de_sat_dt_ice_sclr(ptak) 553 692 !!--------------------------------------------------------------------------------- … … 567 706 de_sat_dt_ice_sclr = LOG(10._wp) * zde * e_sat_ice_sclr(zta) 568 707 END FUNCTION de_sat_dt_ice_sclr 569 !! 708 570 709 FUNCTION de_sat_dt_ice_vctr(ptak) 571 710 !! Same as "e_sat" but over ice rather than water! … … 575 714 !!---------------------------------------------------------------------------------- 576 715 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 577 de_sat_dt_ice_vctr(ji,jj) = de_sat_dt_ice_sclr( ptak(ji,jj) ) 578 END_2D 716 de_sat_dt_ice_vctr(ji,jj) = de_sat_dt_ice_sclr( ptak(ji,jj) ) 717 END_2D 718 579 719 END FUNCTION de_sat_dt_ice_vctr 580 720 581 721 582 583 !===============================================================================================584 585 !===============================================================================================586 722 FUNCTION q_sat_sclr( pta, ppa, l_ice ) 587 723 !!--------------------------------------------------------------------------------- … … 607 743 END IF 608 744 q_sat_sclr = reps0*ze_s/(ppa - (1._wp - reps0)*ze_s) 745 609 746 END FUNCTION q_sat_sclr 610 !! 747 611 748 FUNCTION q_sat_vctr( pta, ppa, l_ice ) 749 612 750 REAL(wp), DIMENSION(jpi,jpj) :: q_sat_vctr 613 751 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute temperature of air [K] … … 620 758 IF( PRESENT(l_ice) ) lice = l_ice 621 759 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 622 q_sat_vctr(ji,jj) = q_sat_sclr( pta(ji,jj) , ppa(ji,jj), l_ice=lice ) 623 END_2D 760 q_sat_vctr(ji,jj) = q_sat_sclr( pta(ji,jj) , ppa(ji,jj), l_ice=lice ) 761 END_2D 762 624 763 END FUNCTION q_sat_vctr 625 !=============================================================================================== 626 627 628 !=============================================================================================== 764 765 629 766 FUNCTION dq_sat_dt_ice_sclr( pta, ppa ) 630 767 !!--------------------------------------------------------------------------------- … … 647 784 ! 648 785 END FUNCTION dq_sat_dt_ice_sclr 649 !! 786 650 787 FUNCTION dq_sat_dt_ice_vctr( pta, ppa ) 788 651 789 REAL(wp), DIMENSION(jpi,jpj) :: dq_sat_dt_ice_vctr 652 790 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute temperature of air [K] … … 655 793 !!---------------------------------------------------------------------------------- 656 794 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 657 dq_sat_dt_ice_vctr(ji,jj) = dq_sat_dt_ice_sclr( pta(ji,jj) , ppa(ji,jj) ) 658 END_2D 795 dq_sat_dt_ice_vctr(ji,jj) = dq_sat_dt_ice_sclr( pta(ji,jj) , ppa(ji,jj) ) 796 END_2D 797 659 798 END FUNCTION dq_sat_dt_ice_vctr 660 !=============================================================================================== 661 662 663 !=============================================================================================== 799 800 664 801 FUNCTION q_air_rh(prha, ptak, ppa) 665 802 !!---------------------------------------------------------------------------------- … … 678 815 ! 679 816 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 680 ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))681 q_air_rh(ji,jj) = ze*reps0/(ppa(ji,jj) - (1. - reps0)*ze)817 ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj)) 818 q_air_rh(ji,jj) = ze*reps0/(ppa(ji,jj) - (1. - reps0)*ze) 682 819 END_2D 683 820 ! … … 685 822 686 823 687 SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, ppa, prlw, &824 SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, ppa, prlw, prhoa, & 688 825 & pQns, pTau, & 689 826 & Qlat) … … 705 842 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa] 706 843 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prlw ! downwelling longwave radiative flux [W/m^2] 844 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prhoa ! air density [kg/m3] 707 845 ! 708 846 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]] … … 715 853 !!---------------------------------------------------------------------------------- 716 854 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 717 zdt = pTa(ji,jj) - pTs(ji,jj) ; zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt ) 718 zdq = pqa(ji,jj) - pqs(ji,jj) ; zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq ) 719 zz0 = pust(ji,jj)/pUb(ji,jj) 720 zCd = zz0*zz0 721 zCh = zz0*ptst(ji,jj)/zdt 722 zCe = zz0*pqst(ji,jj)/zdq 723 724 CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, & 725 & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), & 726 & pTau(ji,jj), zQsen, zQlat ) 727 728 zQlw = qlw_net_sclr( prlw(ji,jj), pTs(ji,jj) ) ! Net longwave flux 729 730 pQns(ji,jj) = zQlat + zQsen + zQlw 731 732 IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 733 END_2D 855 856 zdt = pTa(ji,jj) - pTs(ji,jj) ; zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt ) 857 zdq = pqa(ji,jj) - pqs(ji,jj) ; zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq ) 858 zz0 = pust(ji,jj)/pUb(ji,jj) 859 zCd = zz0*zz0 860 zCh = zz0*ptst(ji,jj)/zdt 861 zCe = zz0*pqst(ji,jj)/zdq 862 863 CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, & 864 & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj), & 865 & pTau(ji,jj), zQsen, zQlat ) 866 867 zQlw = qlw_net_sclr( prlw(ji,jj), pTs(ji,jj) ) ! Net longwave flux 868 869 pQns(ji,jj) = zQlat + zQsen + zQlw 870 871 IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 872 873 END_2D 874 734 875 END SUBROUTINE UPDATE_QNSOL_TAU 735 876 … … 737 878 SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, & 738 879 & pCd, pCh, pCe, & 739 & pwnd, pUb, ppa, 880 & pwnd, pUb, ppa, prhoa, & 740 881 & pTau, pQsen, pQlat, & 741 & pEvap, p rhoa, pfact_evap)742 !!---------------------------------------------------------------------------------- 743 REAL(wp), INTENT(in) :: pzu! height above the sea-level where all this takes place (normally 10m)744 REAL(wp), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]745 REAL(wp), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]746 REAL(wp), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]747 REAL(wp), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]882 & pEvap, pfact_evap ) 883 !!---------------------------------------------------------------------------------- 884 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m) 885 REAL(wp), INTENT(in) :: pTs ! water temperature at the air-sea interface [K] 886 REAL(wp), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg] 887 REAL(wp), INTENT(in) :: pTa ! potential air temperature at z=pzu [K] 888 REAL(wp), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg] 748 889 REAL(wp), INTENT(in) :: pCd 749 890 REAL(wp), INTENT(in) :: pCh 750 891 REAL(wp), INTENT(in) :: pCe 751 REAL(wp), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s] 752 REAL(wp), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 753 REAL(wp), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa] 892 REAL(wp), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s] 893 REAL(wp), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 894 REAL(wp), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa] 895 REAL(wp), INTENT(in) :: prhoa ! Air density at z=pzu [kg/m^3] 754 896 !! 755 897 REAL(wp), INTENT(out) :: pTau ! module of the wind stress [N/m^2] … … 758 900 !! 759 901 REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 760 REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]761 902 REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent) 762 903 !! … … 767 908 IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap 768 909 769 !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 770 ztaa = pTa ! first guess... 771 DO jq = 1, 4 772 zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa ) !#LB: why not "0.5*(pqs+pqa)" rather then "pqa" ??? 773 ztaa = pTa - zgamma*pzu ! Absolute temp. is slightly colder... 774 END DO 775 zrho = rho_air(ztaa, pqa, ppa) 776 zrho = rho_air(ztaa, pqa, ppa-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 777 778 zUrho = pUb*MAX(zrho, 1._wp) ! rho*U10 910 zUrho = pUb*MAX(prhoa, 1._wp) ! rho*U10 779 911 780 912 pTau = zUrho * pCd * pwnd ! Wind stress module … … 785 917 786 918 IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap 787 IF( PRESENT(prhoa) ) prhoa = zrho788 919 789 920 END SUBROUTINE BULK_FORMULA_SCLR 790 !! 921 791 922 SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, & 792 923 & pCd, pCh, pCe, & 793 & pwnd, pUb, ppa, 924 & pwnd, pUb, ppa, prhoa, & 794 925 & pTau, pQsen, pQlat, & 795 & pEvap, p rhoa, pfact_evap )796 !!---------------------------------------------------------------------------------- 797 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)798 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]799 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]800 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]801 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]926 & pEvap, pfact_evap ) 927 !!---------------------------------------------------------------------------------- 928 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m) 929 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTs ! water temperature at the air-sea interface [K] 930 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg] 931 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! potential air temperature at z=pzu [K] 932 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg] 802 933 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd 803 934 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCh 804 935 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCe 805 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s] 806 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 807 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa] 936 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s] 937 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 938 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa] 939 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prhoa ! Air density at z=pzu [kg/m^3] 808 940 !! 809 941 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2] … … 812 944 !! 813 945 REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 814 REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]815 946 REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent) 816 947 !! … … 823 954 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 824 955 825 CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &826 & pCd(ji,jj), pCh(ji,jj), pCe(ji,jj), &827 & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj),&828 & pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj), &829 & pEvap=zevap, prhoa=zrho, pfact_evap=zfact_evap)830 831 IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap832 IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho833 END_2D 956 CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), & 957 & pCd(ji,jj), pCh(ji,jj), pCe(ji,jj), & 958 & pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj), & 959 & pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj), & 960 & pEvap=zevap, pfact_evap=zfact_evap ) 961 962 IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap 963 END_2D 964 834 965 END SUBROUTINE BULK_FORMULA_VCTR 835 966 … … 847 978 !!---------------------------------------------------------------------------------- 848 979 alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79 980 849 981 END FUNCTION alpha_sw_vctr 850 982 … … 861 993 !!---------------------------------------------------------------------------------- 862 994 alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79 995 863 996 END FUNCTION alpha_sw_sclr 864 997 865 998 866 !===============================================================================================867 999 FUNCTION qlw_net_sclr( pdwlw, pts, l_ice ) 868 1000 !!--------------------------------------------------------------------------------- … … 887 1019 zt2 = pts*pts 888 1020 qlw_net_sclr = zemiss*( pdwlw - stefan*zt2*zt2) ! zemiss used both as the IR albedo and IR emissivity... 1021 889 1022 END FUNCTION qlw_net_sclr 890 !! 1023 891 1024 FUNCTION qlw_net_vctr( pdwlw, pts, l_ice ) 1025 892 1026 REAL(wp), DIMENSION(jpi,jpj) :: qlw_net_vctr 893 1027 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdwlw !: downwelling longwave (aka infrared, aka thermal) radiation [W/m^2] … … 900 1034 IF( PRESENT(l_ice) ) lice = l_ice 901 1035 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 902 qlw_net_vctr(ji,jj) = qlw_net_sclr( pdwlw(ji,jj) , pts(ji,jj), l_ice=lice ) 903 END_2D 1036 qlw_net_vctr(ji,jj) = qlw_net_sclr( pdwlw(ji,jj) , pts(ji,jj), l_ice=lice ) 1037 END_2D 1038 904 1039 END FUNCTION qlw_net_vctr 905 !=============================================================================================== 1040 906 1041 907 1042 FUNCTION z0_from_Cd( pzu, pCd, ppsi ) 1043 908 1044 REAL(wp), DIMENSION(jpi,jpj) :: z0_from_Cd !: roughness length [m] 909 1045 REAL(wp) , INTENT(in) :: pzu !: reference height zu [m] … … 921 1057 z0_from_Cd = pzu * EXP( - vkarmn/SQRT(pCd(:,:)) ) !LB: ok, double-checked! 922 1058 END IF 1059 923 1060 END FUNCTION z0_from_Cd 924 1061 1062 925 1063 FUNCTION Cd_from_z0( pzu, pz0, ppsi ) 1064 926 1065 REAL(wp), DIMENSION(jpi,jpj) :: Cd_from_z0 !: (neutral or non-neutral) drag coefficient [] 927 1066 REAL(wp) , INTENT(in) :: pzu !: reference height zu [m] … … 940 1079 END IF 941 1080 Cd_from_z0 = vkarmn2 * Cd_from_z0 * Cd_from_z0 1081 942 1082 END FUNCTION Cd_from_z0 943 1083 … … 965 1105 ! 966 1106 END FUNCTION f_m_louis_sclr 967 !! 1107 968 1108 FUNCTION f_m_louis_vctr( pzu, pRib, pCdn, pz0 ) 1109 969 1110 REAL(wp), DIMENSION(jpi,jpj) :: f_m_louis_vctr 970 1111 REAL(wp), INTENT(in) :: pzu … … 973 1114 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0 974 1115 INTEGER :: ji, jj 975 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 976 f_m_louis_vctr(ji,jj) = f_m_louis_sclr( pzu, pRib(ji,jj), pCdn(ji,jj), pz0(ji,jj) ) 977 END_2D 1116 1117 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1118 f_m_louis_vctr(ji,jj) = f_m_louis_sclr( pzu, pRib(ji,jj), pCdn(ji,jj), pz0(ji,jj) ) 1119 END_2D 1120 978 1121 END FUNCTION f_m_louis_vctr 979 1122 … … 1001 1144 ! 1002 1145 END FUNCTION f_h_louis_sclr 1003 !! 1146 1004 1147 FUNCTION f_h_louis_vctr( pzu, pRib, pChn, pz0 ) 1148 1005 1149 REAL(wp), DIMENSION(jpi,jpj) :: f_h_louis_vctr 1006 1150 REAL(wp), INTENT(in) :: pzu … … 1009 1153 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0 1010 1154 INTEGER :: ji, jj 1011 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1012 f_h_louis_vctr(ji,jj) = f_h_louis_sclr( pzu, pRib(ji,jj), pChn(ji,jj), pz0(ji,jj) ) 1013 END_2D 1155 1156 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1157 f_h_louis_vctr(ji,jj) = f_h_louis_sclr( pzu, pRib(ji,jj), pChn(ji,jj), pz0(ji,jj) ) 1158 END_2D 1159 1014 1160 END FUNCTION f_h_louis_vctr 1161 1015 1162 1016 1163 FUNCTION UN10_from_ustar( pzu, pUzu, pus, ppsi ) … … 1123 1270 1124 1271 1125 !!====================================================================== 1272 1126 1273 END MODULE sbc_phy
Note: See TracChangeset
for help on using the changeset viewer.
