Changeset 11804


Ignore:
Timestamp:
2019-10-25T17:25:19+02:00 (11 months ago)
Author:
laurent
Message:

Use of "TURB_FLUXES@…" inside sbcblk.F90, positive latent HF & positive cool-skin temperature now allowed!

Location:
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90

    r11772 r11804  
    520520      !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 
    521521      !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
    522       !!    (since reanalysis products provide T at z, not theta !) 
     522      !!    (since reanalysis products provide absolute temperature "T" at z, not theta !) 
    523523      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt 
    524524 
     
    594594      END IF ! IF ( ln_skin_cs .OR. ln_skin_wl ) 
    595595 
    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 IF 
    605  
    606596      !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
    607597      !!      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 
    610620         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) 
    613622            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
    614623            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
    615624         END DO 
    616625      END DO 
    617  
    618626      !                          ! add the HF tau contribution to the wind stress module 
    619627      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 
     
    634642      CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
    635643 
    636       !  Turbulent fluxes over ocean 
    637       ! ----------------------------- 
    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 ) THEN 
    643          !! 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 speed 
    645          zqsb (:,:) = cp_air(zqair(:,:)) * zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)  ) ! Sensible Heat, using bulk wind speed 
    646       ELSE 
    647          !! 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 speed 
    650          zqsb (:,:) = cp_air(zqair(:,:)) * zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
    651       ENDIF 
    652  
    653       zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux 
    654  
    655644 
    656645      ! ----------------------------------------------------------------------------- ! 
     
    658647      ! ----------------------------------------------------------------------------- ! 
    659648 
    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.) 
    661650      zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux 
    662651 
     
    681670         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    682671      ! 
    683       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
     672      qns(:,:) = zqlw(:,:) + zqsb(:,:) + zqla(:,:)                                &   ! Downward Non Solar 
    684673         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
    685674         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST !LB??? pst is Celsius !? 
     
    691680      ! 
    692681#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) 
    694683      qsr_oce(:,:) = qsr(:,:) 
    695684#endif 
     
    698687      CALL iom_put( "rho_air"  ,   rhoa )                 ! output air density (kg/m^3) !#LB 
    699688      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 ocean 
    701       CALL iom_put( "qla_oce"  , - zqla )                 ! output downward latent   heat over the ocean 
     689      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 
    702691      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 ocean 
     692      CALL iom_put( "qemp_oce" ,   qns-zqlw-zqsb-zqla )   ! output downward heat content of E-P over the ocean 
    704693      CALL iom_put( "qns_oce"  ,   qns  )                 ! output downward non solar heat over the ocean 
    705694      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  
    495495      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat 
    496496      ! 
    497       REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zUrho, zTs2, zz0, & 
     497      REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zTs2, zz0, & 
    498498         &        zQlat, zQsen, zQlw 
    499499      INTEGER  ::   ji, jj     ! dummy loop indices 
     
    508508            zCh = zz0*ptst(ji,jj)/zdt 
    509509            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             
    514511            CALL TURB_FLUXES( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, & 
    515512               &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj), & 
    516513               &              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)             
    525516            zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux 
    526517 
     
    533524 
    534525 
    535  
    536  
    537  
    538526   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 ) 
    540528      !!---------------------------------------------------------------------------------- 
    541529      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m) 
     
    555543      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat !  [W/m^2] 
    556544      !! 
    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] 
    558547      !! 
    559548      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap 
     
    576565            pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module 
    577566 
    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) =      zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj)) 
    580             pQlat(ji,jj) =  L_vap(pTs(ji,jj)) * zevap 
     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 
    581570 
    582571            IF ( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 
     572            IF ( PRESENT(prhoa) ) prhoa(ji,jj) = zrho 
    583573 
    584574         END DO 
     
    588578 
    589579   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 ) 
    591581      !!---------------------------------------------------------------------------------- 
    592582      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m) 
     
    606596      REAL(wp), INTENT(out) :: pQlat !  [W/m^2] 
    607597      !! 
    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] 
    609600      !! 
    610601      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap 
     
    625616      pTau = zUrho * pCd * pwnd ! Wind stress module 
    626617 
    627       zevap        = MIN( zUrho * pCe * (pqa - pqs) , 0._wp )   ! we do not want condensation & Qlat > 0 ! 
    628       pQsen =      zUrho * pCh * (pTa - pTs) * cp_air(pqa) 
    629       pQlat =  L_vap(pTs) * zevap 
     618      zevap = zUrho * pCe * (pqa - pqs) 
     619      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa) 
     620      pQlat = L_vap(pTs) * zevap 
    630621 
    631622      IF ( PRESENT(pEvap) ) pEvap = - zevap 
     623      IF ( PRESENT(prhoa) ) prhoa = zrho 
    632624 
    633625   END SUBROUTINE TURB_FLUXES_SCLR 
    634  
    635626 
    636627 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90

    r11775 r11804  
    6363      !!--------------------------------------------------------------------- 
    6464      !! 
    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) 
    6666      !! 
    6767      !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 
     
    7878      !!     *pSST*       bulk SST (taken at depth gdept_1d(1))          [K] 
    7979      !!     *pQlat*      surface latent heat flux                       [K] 
    80       !! 
    8180      !!------------------------------------------------------------------ 
    82       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQsw   ! net solar a.k.a shortwave radiation into the ocean (after albedo) [W/m^2] 
    83       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQnsol ! non-solar heat flux to the ocean [W/m^2] 
    84       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pustar  ! friction velocity, temperature and humidity (u*,t*,q*) 
    85       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pSST ! bulk SST [K] 
    86       REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pQlat  ! latent heat flux [W/m^2] 
     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] 
    8786      !!--------------------------------------------------------------------- 
    8887      INTEGER  :: ji, jj, jc 
     
    9291         DO ji = 1, jpi 
    9392 
    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... 
    9695 
    9796            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 
    9897 
    9998            DO jc = 1, 4 ! because implicit in terms of zdelta... 
    100                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 ! 
    101                zQabs  = MIN( -0.1_wp , pQnsol(ji,jj) + zfr*pQsw(ji,jj) ) ! Total cooling at the interface 
     99               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) 
    102101               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 
    103102            END DO 
    104103 
    105             dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp )   ! temperature increment 
     104            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 
    106105 
    107106         END DO 
     
    189188                  &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 
    190189               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=', zfr 
    193190 
    194191               IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN 
     
    196193                  ! since zQabs <= 0._wp 
    197194                  ! => 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=', zQabs 
    200                   !PRINT *, '#LBD: => leaving without changing anything...' 
    201195                  l_exit = .TRUE. 
    202196               END IF 
     
    207201            !LOLO: remove??? has a strong influence !!! 
    208202            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*rdt 
    210                !PRINT *, '#LBD:  => time to destroy the warm-layer!' 
    211203               l_exit       = .TRUE. 
    212204               l_destroy_wl = .TRUE. 
     
    219211               ! 1/ A warm layer already exists (dT>0) but it is cooling down because Qabs<0 
    220212               ! 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) 
    226213 
    227214               ztac = Tau_ac(ji,jj) + MAX(.002_wp , pTau(ji,jj))*rdt      ! updated momentum integral 
     
    238225                  zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth 
    239226               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' 
    243227 
    244228               IF ( zqac <= 0._wp ) THEN 
     
    263247            END IF 
    264248 
    265             !PRINT *, '#LBD: exit values for Qac & Tac:', REAL(zqac,4), REAL(ztac,4) 
    266  
    267249            IF ( iwait == 0 ) THEN 
    268250               !! Iteration loop within bulk algo is over, time to update what needs to be updated: 
    269251               dT_wl(ji,jj)  = zdTwl 
    270252               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) 
    272253               Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral 
    273254               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' 
    276255            END IF 
    277256 
     
    284263 
    285264 
    286    FUNCTION delta_skin_layer( palpha, pQabs, pQlat, pustar_a ) 
     265   FUNCTION delta_skin_layer( palpha, pQd, pQlat, pustar_a ) 
    287266      !!--------------------------------------------------------------------- 
    288267      !! Computes the thickness (m) of the viscous skin layer. 
     
    298277      REAL(wp)                :: delta_skin_layer 
    299278      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!) 
    300       REAL(wp), INTENT(in)    :: pQabs    ! < 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 
     279      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 
    301280      REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2] 
    302281      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s] 
    303282      !!--------------------------------------------------------------------- 
    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      ! 
    309291      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water 
    310292      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 
    316300   END FUNCTION delta_skin_layer 
    317301 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90

    r11772 r11804  
    8888         DO ji = 1, jpi 
    8989 
    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... 
    9392 
    9493            zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 
    9594 
    9695            DO jc = 1, 4 ! because implicit in terms of zdelta... 
    97                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 
     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 
    9897               !              =>  (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) 
    10199               zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 
    102100            END DO 
    103101 
    104             dT_cs(ji,jj) = MIN( zQabs*zdelta/rk0_w , 0._wp )   ! temperature increment 
     102            dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 
    105103 
    106104         END DO 
     
    142140         & zalpha, & !: thermal expansion coefficient of sea-water [1/K] 
    143141         & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL 
    144          & zfr, zeta, ztmp, & 
     142         & zfr, zeta, & 
    145143         & zusw, zusw2, & 
    146144         & zLa, zfLa, & 
     
    239237 
    240238 
    241    FUNCTION delta_skin_layer( palpha, pQabs, pustar_a ) 
     239   FUNCTION delta_skin_layer( palpha, pQd, pustar_a ) 
    242240      !!--------------------------------------------------------------------- 
    243241      !! Computes the thickness (m) of the viscous skin layer. 
     
    253251      REAL(wp)                :: delta_skin_layer 
    254252      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!) 
    255       REAL(wp), INTENT(in)    :: pQabs    ! < 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 
     253      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 
    256254      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s] 
    257255      !!--------------------------------------------------------------------- 
    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      ! 
    267261      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water 
    268262      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 
    277270   END FUNCTION delta_skin_layer 
    278  
    279271 
    280272 
Note: See TracChangeset for help on using the changeset viewer.