Changeset 11804 for NEMO/branches
- Timestamp:
- 2019-10-25T17:25:19+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC
- Files:
-
- 4 edited
-
sbcblk.F90 (modified) (7 diffs)
-
sbcblk_phy.F90 (modified) (8 diffs)
-
sbcblk_skin_coare.F90 (modified) (11 diffs)
-
sbcblk_skin_ecmwf.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90
r11772 r11804 520 520 !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 521 521 !! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 522 !! (since reanalysis products provide Tat z, not theta !)522 !! (since reanalysis products provide absolute temperature "T" at z, not theta !) 523 523 ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt 524 524 … … 594 594 END IF ! IF ( ln_skin_cs .OR. ln_skin_wl ) 595 595 596 597 598 599 ! ! Compute true air density :600 IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove rhoa*grav*rn_zu from SLP...)601 rhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) )602 ELSE ! At zt:603 rhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), zqair(:,:), sf(jp_slp)%fnow(:,:,1) )604 END IF605 606 596 !! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef. 607 597 !! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef. 608 609 DO jj = 1, jpj ! tau module, i and j component 598 599 IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 600 !! If zu == zt, then ensuring once for all that: 601 t_zu(:,:) = ztpot(:,:) 602 q_zu(:,:) = zqair(:,:) 603 END IF 604 605 606 ! Turbulent fluxes over ocean => TURB_FLUXES @ sbcblk_phy.F90 607 ! ------------------------------------------------------------- 608 609 CALL TURB_FLUXES( rn_zu, zst(:,:), zsq(:,:), t_zu(:,:), q_zu(:,:), Cd_atm(:,:), Ch_atm(:,:), Ce_atm(:,:), & 610 & wndm(:,:), zU_zu(:,:), sf(jp_slp)%fnow(:,:,1), & 611 & taum(:,:), zqsb(:,:), zqla(:,:), & 612 & pEvap=zevap(:,:), prhoa=rhoa(:,:) ) 613 614 zqla(:,:) = zqla(:,:) * tmask(:,:,1) 615 zqsb(:,:) = zqsb(:,:) * tmask(:,:,1) 616 taum(:,:) = taum(:,:) * tmask(:,:,1) 617 zevap(:,:) = zevap(:,:) * tmask(:,:,1) 618 619 DO jj = 1, jpj ! tau i and j component on T-grid points 610 620 DO ji = 1, jpi 611 zztmp = rhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed 612 taum (ji,jj) = zztmp * wndm (ji,jj) 621 zztmp = taum(ji,jj) / wndm(ji,jj) 613 622 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 614 623 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 615 624 END DO 616 625 END DO 617 618 626 ! ! add the HF tau contribution to the wind stress module 619 627 IF( lhftau ) taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) … … 634 642 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 635 643 636 ! Turbulent fluxes over ocean637 ! -----------------------------638 639 ! zqla used as temporary array, for rho*U (common term of bulk formulae):640 zqla(:,:) = rhoa(:,:) * zU_zu(:,:) * tmask(:,:,1)641 642 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN643 !! q_air and t_air are given at 10m (wind reference height)644 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - zqair(:,:)) ) ! Evaporation, using bulk wind speed645 zqsb (:,:) = cp_air(zqair(:,:)) * zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed646 ELSE647 !! q_air and t_air are not given at 10m (wind reference height)648 ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!649 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:)) ) ! Evaporation, using bulk wind speed650 zqsb (:,:) = cp_air(zqair(:,:)) * zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed651 ENDIF652 653 zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:) ! Latent Heat flux654 655 644 656 645 ! ----------------------------------------------------------------------------- ! … … 658 647 ! ----------------------------------------------------------------------------- ! 659 648 660 !! LB: now after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)649 !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST ! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 661 650 zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux 662 651 … … 681 670 & - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1) 682 671 ! 683 qns(:,:) = zqlw(:,:) - zqsb(:,:) -zqla(:,:) & ! Downward Non Solar672 qns(:,:) = zqlw(:,:) + zqsb(:,:) + zqla(:,:) & ! Downward Non Solar 684 673 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 685 674 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST !LB??? pst is Celsius !? … … 691 680 ! 692 681 #if defined key_si3 693 qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) -zqla(:,:) ! non solar without emp (only needed by SI3)682 qns_oce(:,:) = zqlw(:,:) + zqsb(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3) 694 683 qsr_oce(:,:) = qsr(:,:) 695 684 #endif … … 698 687 CALL iom_put( "rho_air" , rhoa ) ! output air density (kg/m^3) !#LB 699 688 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 700 CALL iom_put( "qsb_oce" , -zqsb ) ! output downward sensible heat over the ocean701 CALL iom_put( "qla_oce" , -zqla ) ! output downward latent heat over the ocean689 CALL iom_put( "qsb_oce" , zqsb ) ! output downward sensible heat over the ocean 690 CALL iom_put( "qla_oce" , zqla ) ! output downward latent heat over the ocean 702 691 CALL iom_put( "evap_oce" , zevap ) ! evaporation 703 CALL iom_put( "qemp_oce" , qns-zqlw +zqsb+zqla ) ! output downward heat content of E-P over the ocean692 CALL iom_put( "qemp_oce" , qns-zqlw-zqsb-zqla ) ! output downward heat content of E-P over the ocean 704 693 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 705 694 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90
r11772 r11804 495 495 REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat 496 496 ! 497 REAL(wp) :: zdt, zdq, zCd, zCh, zCe, z Urho, zTs2, zz0, &497 REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zTs2, zz0, & 498 498 & zQlat, zQsen, zQlw 499 499 INTEGER :: ji, jj ! dummy loop indices … … 508 508 zCh = zz0*ptst(ji,jj)/zdt 509 509 zCe = zz0*pqst(ji,jj)/zdq 510 511 !zUrho = pUb(ji,jj)*MAX(rho_air(pTa(ji,jj), pqa(ji,jj), pslp(ji,jj)), 1._wp) ! rho*U10 512 zTs2 = pTs(ji,jj)*pTs(ji,jj) 513 510 514 511 CALL TURB_FLUXES( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, & 515 512 & pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj), & 516 513 & pTau(ji,jj), zQsen, zQlat ) 517 518 519 ! Wind stress module: 520 !pTau(ji,jj) = zCd*zUrho*pUb(ji,jj) ! lolo? 521 522 ! Non-Solar heat flux to the ocean: 523 !zQlat = MIN ( zUrho*zCe*L_vap( pTs(ji,jj)) * zdq , 0._wp ) ! we do not want a Qlat > 0 ! 524 !zQsen = zUrho*zCh*cp_air(pqa(ji,jj)) * zdt 514 515 zTs2 = pTs(ji,jj)*pTs(ji,jj) 525 516 zQlw = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux 526 517 … … 533 524 534 525 535 536 537 538 526 SUBROUTINE TURB_FLUXES_VCTR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, & 539 & pTau, pQsen, pQlat, pEvap )527 & pTau, pQsen, pQlat, pEvap, prhoa ) 540 528 !!---------------------------------------------------------------------------------- 541 529 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m) … … 555 543 REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat ! [W/m^2] 556 544 !! 557 REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! [kg/m^2/s] 545 REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 546 REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3] 558 547 !! 559 548 REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap … … 576 565 pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module 577 566 578 zevap = MIN( zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj)) , 0._wp ) ! we do not want condensation & Qlat > 0 !579 pQsen(ji,jj) = 580 pQlat(ji,jj) = 567 zevap = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj)) 568 pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj)) 569 pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap 581 570 582 571 IF ( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 572 IF ( PRESENT(prhoa) ) prhoa(ji,jj) = zrho 583 573 584 574 END DO … … 588 578 589 579 SUBROUTINE TURB_FLUXES_SCLR( pzu, pTs, pqs, pTa, pqa, pCd, pCh, pCe, pwnd, pUb, pslp, & 590 & pTau, pQsen, pQlat, pEvap )580 & pTau, pQsen, pQlat, pEvap, prhoa ) 591 581 !!---------------------------------------------------------------------------------- 592 582 REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m) … … 606 596 REAL(wp), INTENT(out) :: pQlat ! [W/m^2] 607 597 !! 608 REAL(wp), INTENT(out), OPTIONAL :: pEvap ! [kg/m^2/s] 598 REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 599 REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3] 609 600 !! 610 601 REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap … … 625 616 pTau = zUrho * pCd * pwnd ! Wind stress module 626 617 627 zevap = MIN( zUrho * pCe * (pqa - pqs) , 0._wp ) ! we do not want condensation & Qlat > 0 !628 pQsen = 629 pQlat = 618 zevap = zUrho * pCe * (pqa - pqs) 619 pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa) 620 pQlat = L_vap(pTs) * zevap 630 621 631 622 IF ( PRESENT(pEvap) ) pEvap = - zevap 623 IF ( PRESENT(prhoa) ) prhoa = zrho 632 624 633 625 END SUBROUTINE TURB_FLUXES_SCLR 634 635 626 636 627 -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90
r11775 r11804 63 63 !!--------------------------------------------------------------------- 64 64 !! 65 !! Cool-Skin scheme according to Fairall et al.1996, revisited for COARE 3.6 (Fairall et al., 2019)65 !! Cool-skin parameterization, based on Fairall et al., 1996, revisited for COARE 3.6 (Fairall et al., 2019) 66 66 !! 67 67 !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., … … 78 78 !! *pSST* bulk SST (taken at depth gdept_1d(1)) [K] 79 79 !! *pQlat* surface latent heat flux [K] 80 !!81 80 !!------------------------------------------------------------------ 82 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) 83 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) 84 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) 85 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) 86 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) 81 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQsw ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 82 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 83 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pustar ! friction velocity, temperature and humidity (u*,t*,q*) 84 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pSST ! bulk SST [K] 85 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pQlat ! latent heat flux [W/m^2] 87 86 !!--------------------------------------------------------------------- 88 87 INTEGER :: ji, jj, jc … … 92 91 DO ji = 1, jpi 93 92 94 zQabs = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$95 ! ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$93 zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta, 94 ! ! => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta... 96 95 97 96 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 98 97 99 98 DO jc = 1, 4 ! because implicit in terms of zdelta... 100 zfr 101 zQabs = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface99 zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) / !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 ! 100 zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 102 101 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 103 102 END DO 104 103 105 dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp ) ! temperature increment104 dT_cs(ji,jj) = zQabs*zdelta/rk0_w ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 106 105 107 106 END DO … … 189 188 & + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 190 189 zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer !LOLO: depends of zfr, which is wild guess... Wrong!!! 191 !PRINT *, '#LBD: Initial Qsw & Qnsol:', NINT(pQsw(ji,jj)), NINT(pQnsol(ji,jj))192 !PRINT *, '#LBD: =>Qabs:', zQabs,' zfr=', zfr193 190 194 191 IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN … … 196 193 ! since zQabs <= 0._wp 197 194 ! => no need to go further 198 !PRINT *, '#LBD: we have not started to to build a WL yet (dT==0)'199 !PRINT *, '#LBD: and theres no way it can occur now since zQabs=', zQabs200 !PRINT *, '#LBD: => leaving without changing anything...'201 195 l_exit = .TRUE. 202 196 END IF … … 207 201 !LOLO: remove??? has a strong influence !!! 208 202 IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN 209 !PRINT *, '#LBD: Oh boy! Next Qnt_ac looking weak! =>', Qnt_ac(ji,jj) + zQabs*rdt210 !PRINT *, '#LBD: => time to destroy the warm-layer!'211 203 l_exit = .TRUE. 212 204 l_destroy_wl = .TRUE. … … 219 211 ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0 220 212 ! 2/ Regardless of WL formed (dT==0 or dT>0), we are in the process to initiate one or warm further it ! 221 222 !PRINT *, '#LBD:======================================================'223 !PRINT *, '#LBD: WL action makes sense now! => zQabs,dT_wl=', REAL(zQabs,4), REAL(zdTwl,4)224 !PRINT *, '#LBD:======================================================'225 !PRINT *, '#LBD: current values for Qac and Tac=', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4)226 213 227 214 ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt ! updated momentum integral … … 238 225 zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth 239 226 END DO 240 !PRINT *, '#LBD: updated absorption and WL depth=', REAL(zfr,4), REAL(zHwl,4)241 !PRINT *, '#LBD: updated value for Qabs=', REAL(zQabs,4), 'W/m2'242 !PRINT *, '#LBD: updated value for Qac =', REAL(zqac,4), 'J'243 227 244 228 IF ( zqac <= 0._wp ) THEN … … 263 247 END IF 264 248 265 !PRINT *, '#LBD: exit values for Qac & Tac:', REAL(zqac,4), REAL(ztac,4)266 267 249 IF ( iwait == 0 ) THEN 268 250 !! Iteration loop within bulk algo is over, time to update what needs to be updated: 269 251 dT_wl(ji,jj) = zdTwl 270 252 Hz_wl(ji,jj) = zHwl 271 !PRINT *, '#LBD: FINAL EXIT values for dT_wl & Hz_wl:', REAL(dT_wl(ji,jj),4), REAL(Hz_wl(ji,jj),4)272 253 Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral 273 254 Tau_ac(ji,jj) = ztac 274 !PRINT *, '#LBD: FINAL EXIT values for Qac & Tac:', REAL(Qnt_ac(ji,jj),4), REAL(Tau_ac(ji,jj),4)275 !PRINT *, '#LBD'276 255 END IF 277 256 … … 284 263 285 264 286 FUNCTION delta_skin_layer( palpha, pQ abs, pQlat, pustar_a )265 FUNCTION delta_skin_layer( palpha, pQd, pQlat, pustar_a ) 287 266 !!--------------------------------------------------------------------- 288 267 !! Computes the thickness (m) of the viscous skin layer. … … 298 277 REAL(wp) :: delta_skin_layer 299 278 REAL(wp), INTENT(in) :: palpha ! thermal expansion coefficient of sea-water (SST accurate enough!) 300 REAL(wp), INTENT(in) :: pQ abs! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996279 REAL(wp), INTENT(in) :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 301 280 REAL(wp), INTENT(in) :: pQlat ! latent heat flux [W/m^2] 302 281 REAL(wp), INTENT(in) :: pustar_a ! friction velocity in the air (u*) [m/s] 303 282 !!--------------------------------------------------------------------- 304 REAL(wp) :: zusw, zusw2, zlamb, zQb 305 !!--------------------------------------------------------------------- 306 307 zQb = pQabs + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha ! LOLO: Double check sign + division by palpha !!! units are okay! 308 283 REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp 284 !!--------------------------------------------------------------------- 285 286 zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha ! LOLO: Double check sign + division by palpha !!! units are okay! 287 288 ztf = 0.5_wp + SIGN(0.5_wp, zQd) ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case) 289 ! ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux) 290 ! 309 291 zusw = MAX(pustar_a, 1.E-4_wp) * sq_radrw ! u* in the water 310 292 zusw2 = zusw*zusw 311 312 zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 313 314 delta_skin_layer = zlamb*rnu0_w/zusw 315 293 ! 294 zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*zQd, 0._wp)**0.75 )**(-1./3.) ! see Eq.(14) in Fairall et al., 1996 295 ! => zlamb is not used when Qd > 0, and since zcon0 < 0, we just use this "MAX" to prevent FPE errors (something_negative)**0.75 296 ! 297 ztmp = rnu0_w/zusw 298 delta_skin_layer = (1._wp-ztf) * zlamb*ztmp & ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996 299 & + ztf * MIN(6._wp*ztmp , 0.007_wp) ! when Qd > 0 316 300 END FUNCTION delta_skin_layer 317 301 -
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90
r11772 r11804 88 88 DO ji = 1, jpi 89 89 90 zQabs = MIN( -0.1_wp , pQnsol(ji,jj) ) ! first guess, we do not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes$ 91 ! ! also, we ONLY consider when the viscous layer is loosing heat to the atmosphere, we only deal with cool-skin! => hence the "MIN( -0$ 92 !zQabs = pQnsol(ji,jj) 90 zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta, 91 ! ! => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta... 93 92 94 93 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 95 94 96 95 DO jc = 1, 4 ! because implicit in terms of zdelta... 97 zfr 96 zfr = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005 98 97 ! => (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996) 99 zQabs = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface 100 !zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 98 zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 101 99 zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 102 100 END DO 103 101 104 dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp ) ! temperature increment102 dT_cs(ji,jj) = zQabs*zdelta/rk0_w ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 105 103 106 104 END DO … … 142 140 & zalpha, & !: thermal expansion coefficient of sea-water [1/K] 143 141 & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL 144 & zfr, zeta, ztmp,&142 & zfr, zeta, & 145 143 & zusw, zusw2, & 146 144 & zLa, zfLa, & … … 239 237 240 238 241 FUNCTION delta_skin_layer( palpha, pQ abs, pustar_a )239 FUNCTION delta_skin_layer( palpha, pQd, pustar_a ) 242 240 !!--------------------------------------------------------------------- 243 241 !! Computes the thickness (m) of the viscous skin layer. … … 253 251 REAL(wp) :: delta_skin_layer 254 252 REAL(wp), INTENT(in) :: palpha ! thermal expansion coefficient of sea-water (SST accurate enough!) 255 REAL(wp), INTENT(in) :: pQ abs! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996253 REAL(wp), INTENT(in) :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 256 254 REAL(wp), INTENT(in) :: pustar_a ! friction velocity in the air (u*) [m/s] 257 255 !!--------------------------------------------------------------------- 258 REAL(wp) :: zusw, zusw2, zlamb, zQb 259 !!--------------------------------------------------------------------- 260 261 zQb = pQabs 262 263 !zQ = MIN( -0.1_wp , pQabs ) 264 265 !ztf = 0.5_wp + SIGN(0.5_wp, zQ) ! Qabs < 0 => cooling of the layer => ztf = 0 (normal case) 266 ! ! Qabs > 0 => warming of the layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux) 256 REAL(wp) :: zusw, zusw2, zlamb, ztf, ztmp 257 !!--------------------------------------------------------------------- 258 ztf = 0.5_wp + SIGN(0.5_wp, pQd) ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case) 259 ! ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux) 260 ! 267 261 zusw = MAX(pustar_a, 1.E-4_wp) * sq_radrw ! u* in the water 268 262 zusw2 = zusw*zusw 269 270 zlamb = 6._wp*( 1._wp + (palpha*zcon0/(zusw2*zusw2)*zQb)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 271 !zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*zQ, 0._wp)**0.75 )**(-1./3.) ! see eq.(14) in Fairall et al., 1996 272 273 delta_skin_layer = zlamb*rnu0_w/zusw 274 275 !delta_skin_layer = (1._wp - ztf) * zlamb*rnu0_w/zusw & ! see eq.(12) in Fairall et al., 1996 276 ! & + ztf * MIN(6._wp*rnu0_w/zusw , 0.007_wp) 263 ! 264 zlamb = 6._wp*( 1._wp + MAX(palpha*zcon0/(zusw2*zusw2)*pQd, 0._wp)**0.75 )**(-1./3.) ! see Eq.(14) in Fairall et al., 1996 265 ! => zlamb is not used when Qd > 0, and since zcon0 < 0, we just use this "MAX" to prevent FPE errors (something_negative)**0.75 266 ! 267 ztmp = rnu0_w/zusw 268 delta_skin_layer = (1._wp-ztf) * zlamb*ztmp & ! regular case, Qd < 0, see Eq.(12) in Fairall et al., 1996 269 & + ztf * MIN(6._wp*ztmp , 0.007_wp) ! when Qd > 0 277 270 END FUNCTION delta_skin_layer 278 279 271 280 272
Note: See TracChangeset
for help on using the changeset viewer.
