Changeset 12313


Ignore:
Timestamp:
2020-01-14T13:40:28+01:00 (2 weeks ago)
Author:
agn
Message:

Modify sbcblk.F90 to accept humidity as dewpoint

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/SBC/sbcblk.F90

    r12178 r12313  
    124124   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 
    125125 
     126   LOGICAL  ::   ln_humi_dpt = .FALSE.                                        ! calculate specific hunidity from dewpoint 
     127   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qair                      ! specific humidity of air at input height 
     128 
    126129   INTEGER  ::   nblk           ! choice of the bulk algorithm 
    127130   !                            ! associated indices: 
     
    145148      !!------------------------------------------------------------------- 
    146149      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 
    147          &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 
     150         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), qair(jpi,jpj), STAT=sbc_blk_alloc ) 
    148151      ! 
    149152      CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) 
     
    171174      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    172175         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                & 
    173          &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
     176         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, ln_humi_dpt,&   ! bulk algorithm 
    174177         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
    175178         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 
     
    323326      ! 
    324327      !                                            ! compute the surface ocean fluxes using bulk formulea 
     328      ! ..... if dewpoint supplied instead of specific humidaity, calculate specific humidity 
     329      IF(ln_humi_dpt) THEN 
     330         qair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     331      ELSE 
     332         qair(:,:) = sf(jp_humi)%fnow(:,:,1) 
     333      END IF 
     334       
    325335      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 
    326336 
     
    332342         ENDIF  
    333343         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    334          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1) 
     344         qatm_ice(:,:)    = qair(:,:) 
    335345         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
    336346         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac 
     
    434444      !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 
    435445      !!    (since reanalysis products provide T at z, not theta !) 
    436       ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt 
     446      ztpot(:,:) = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), qair(:,:) ) * rn_zqt 
    437447 
    438448      SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    439449      ! 
    440       CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2 
     450      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! NCAR-COREv2 
    441451         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    442       CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0 
     452      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! COARE v3.0 
    443453         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    444       CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5 
     454      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! COARE v3.5 
    445455         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    446       CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF 
     456      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! ECMWF 
    447457         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    448458      CASE DEFAULT 
     
    454464         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) ) 
    455465      ELSE                                      ! At zt: 
    456          zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     466         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1) ) 
    457467      END IF 
    458468 
     
    495505      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    496506         !! q_air and t_air are given at 10m (wind reference height) 
    497          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
    498          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
     507         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - qair(:,:)) ) ! Evaporation, using bulk wind speed 
     508         zqsb (:,:) = cp_air(qair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
    499509      ELSE 
    500510         !! q_air and t_air are not given at 10m (wind reference height) 
    501511         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
    502512         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 
    503          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
     513         zqsb (:,:) = cp_air(qair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
    504514      ENDIF 
    505515 
     
    742752      ! local scalars ( place there for vector optimisation purposes) 
    743753      ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    744       zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
     754      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1)) 
    745755 
    746756      !!gm brutal.... 
     
    806816      zcoef_dqla = -Ls * 11637800. * (-5897.8) 
    807817      ! 
    808       zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     818      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1) ) 
    809819      ! 
    810820      zztmp = 1. / ( 1. - albo ) 
     
    837847               ! Latent Heat 
    838848               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
    839                   &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
     849                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - qair(ji,jj) ) ) 
    840850               ! Latent heat sensitivity for ice (Dqla/Dt) 
    841851               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
Note: See TracChangeset for help on using the changeset viewer.