New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12199 for NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/sbcblk.F90 – NEMO

Ignore:
Timestamp:
2019-12-12T09:30:08+01:00 (4 years ago)
Author:
laurent
Message:

Catches up with SBCBLK bugfixes of branch "dev_r12072_MERGE_OPTION2_2019"

File:
1 edited

Legend:

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

    r12193 r12199  
    178178      ! 
    179179      !                             !** read bulk namelist 
     180      REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters 
    180181      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) 
    181182901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' ) 
    182183      ! 
     184      REWIND( numnam_cfg )                !* Namelist namsbc_blk in configuration namelist : bulk parameters 
    183185      READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 ) 
    184186902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' ) 
     
    399401      ! Sanity/consistence test on humidity at first time step to detect potential screw-up: 
    400402      IF( kt == nit000 ) THEN 
    401          WRITE(numout,*) '' 
     403         IF(lwp) WRITE(numout,*) '' 
    402404#if defined key_agrif 
    403          WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ===' 
     405         IF(lwp) WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ===' 
    404406#else 
    405407         ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain 
     
    415417            END SELECT 
    416418            IF(ztmp < 0._wp) THEN 
    417                WRITE(numout,'("   Mean humidity value found on proc #",i5.5," is: ",f)') narea, ztmp 
     419               IF (lwp) WRITE(numout,'("   Mean humidity value found on proc #",i5.5," is: ",f)') narea, ztmp 
    418420               CALL ctl_stop( 'STOP', 'Something is wrong with air humidity!!!', & 
    419421                  &   ' ==> check the unit in your input files'       , & 
     
    422424            END IF 
    423425         END IF 
    424          WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 
     426         IF(lwp) WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 
    425427#endif 
    426          WRITE(numout,*) '' 
     428         IF(lwp) WRITE(numout,*) '' 
    427429      END IF !IF( kt == nit000 ) 
    428430      !                                            ! compute the surface ocean fluxes using bulk formulea 
     
    432434            &                sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m,       &   !   <<= in 
    433435            &                sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1),   &   !   <<= in (wl/cs) 
    434             &                zssq, zcd_du, zsen, zevp )                              !   =>> out 
     436            &                tsk_m, zssq, zcd_du, zsen, zevp )                       !   =>> out 
    435437 
    436438         CALL blk_oce_2(     sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in 
    437439            &                sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1),   &   !   <<= in 
    438             &                sf(jp_snow)%fnow(:,:,1), sst_m,                     &   !   <<= in 
     440            &                sf(jp_snow)%fnow(:,:,1), tsk_m,                     &   !   <<= in 
    439441            &                zsen, zevp )                                            !   <=> in out 
    440442      ENDIF 
     
    470472 
    471473   SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, &  ! inp 
    472       &              pslp , pst   , pu   , pv,    &  ! inp 
    473       &              pqsr , pqlw  ,               &  ! inp 
    474       &              pssq , pcd_du, psen , pevp   )  ! out 
     474      &                  pslp , pst   , pu   , pv,        &  ! inp 
     475      &                  pqsr , pqlw  ,                   &  ! inp 
     476      &                  ptsk, pssq , pcd_du, psen , pevp   )  ! out 
    475477      !!--------------------------------------------------------------------- 
    476478      !!                     ***  ROUTINE blk_oce_1  *** 
     
    494496      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   ptair  ! potential temperature at T-points        [Kelvin] 
    495497      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pslp   ! sea-level pressure                       [Pa] 
    496       REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pst    ! surface temperature                      [Celcius] 
     498      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pst    ! surface temperature                      [Celsius] 
    497499      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s] 
    498500      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s] 
    499501      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqsr   ! 
    500502      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqlw   ! 
     503      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   ptsk   ! skin temp. (or SST if CS & WL not used)  [Celsius] 
    501504      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pssq   ! specific humidity at pst                 [kg/kg] 
    502505      REAL(wp), INTENT(  out), DIMENSION(:,:) ::   pcd_du ! Cd x |dU| at T-points                    [m/s] 
     
    507510      REAL(wp) ::   zztmp                ! local variable 
    508511      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
    509       REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
    510512      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    511513      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
     
    519521      ! 
    520522      ! local scalars ( place there for vector optimisation purposes) 
    521       zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
     523      !                           ! Temporary conversion from Celcius to Kelvin (and set minimum value far above 0 K) 
     524      ptsk(:,:) = pst(:,:) + rt0  ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!) 
    522525 
    523526      ! ----------------------------------------------------------------------------- ! 
     
    566569 
    567570      ! specific humidity at SST 
    568       pssq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), pslp(:,:) ) 
     571      pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) ) 
    569572 
    570573      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    571          zztmp1(:,:) = zst(:,:) 
     574         !! Backup "bulk SST" and associated spec. hum. 
     575         zztmp1(:,:) = ptsk(:,:) 
    572576         zztmp2(:,:) = pssq(:,:) 
    573577      ENDIF 
     
    608612 
    609613      CASE( np_NCAR      ) 
    610          CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm,                              & 
     614         CALL turb_ncar    ( rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm,                              & 
    611615            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    612616 
    613617      CASE( np_COARE_3p0 ) 
    614          CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     618         CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
    615619            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
    616620            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
    617621 
    618622      CASE( np_COARE_3p6 ) 
    619          CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     623         CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
    620624            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
    621625            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
    622626 
    623627      CASE( np_ECMWF     ) 
    624          CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
     628         CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
    625629            &                zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
    626630            &                Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 
     
    632636 
    633637      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    634          !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and pssq: 
    635          WHERE ( fr_i < 0.001_wp ) 
    636             ! zst and pssq have been updated by cool-skin/warm-layer scheme and we keep it!!! 
    637             zst(:,:)  =  zst(:,:)*tmask(:,:,1) 
    638             pssq(:,:) = pssq(:,:)*tmask(:,:,1) 
    639          ELSEWHERE 
    640             ! we forget about the update... 
    641             zst(:,:)  = zztmp1(:,:) !#LB: using what we backed up before skin-algo 
    642             pssq(:,:) = zztmp2(:,:) !#LB:  "   "   " 
     638         !! ptsk and pssq have been updated!!! 
     639         !! 
     640         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq: 
     641         WHERE ( fr_i(:,:) > 0.001_wp ) 
     642            ! sea-ice present, we forget about the update, using what we backed up before call to turb_*() 
     643            ptsk(:,:) = zztmp1(:,:) 
     644            pssq(:,:) = zztmp2(:,:) 
    643645         END WHERE 
    644646      END IF 
     
    669671         END DO 
    670672      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation 
    671          CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 
     673         CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 
    672674            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         & 
    673675            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 & 
     
    707709         ENDIF 
    708710         ! 
    709       ENDIF 
    710       ! 
     711      ENDIF !IF( ln_abl ) 
     712       
     713      ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1)  ! Back to Celsius 
     714             
     715      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
     716         CALL iom_put( "t_skin" ,  ptsk        )  ! T_skin in Celsius 
     717         CALL iom_put( "dt_skin" , ptsk - pst  )  ! T_skin - SST temperature difference... 
     718      ENDIF 
     719 
    711720      IF(ln_ctl) THEN 
    712721         CALL prt_ctl( tab2d_1=pevp  , clinfo1=' blk_oce_1: pevp   : ' ) 
     
    719728 
    720729   SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec,   &   ! <<= in 
    721       &          psnow, pst , psen, pevp     )   ! <<= in 
     730      &                  psnow, ptsk, psen, pevp     )   ! <<= in 
    722731      !!--------------------------------------------------------------------- 
    723732      !!                     ***  ROUTINE blk_oce_2  *** 
     
    740749      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pprec 
    741750      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psnow 
    742       REAL(wp), INTENT(in), DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius] 
     751      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptsk   ! SKIN surface temperature   [Celsius] 
    743752      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psen 
    744753      REAL(wp), INTENT(in), DIMENSION(:,:) ::   pevp 
     
    746755      INTEGER  ::   ji, jj               ! dummy loop indices 
    747756      REAL(wp) ::   zztmp,zz1,zz2,zz3    ! local variable 
    748       REAL(wp), DIMENSION(jpi,jpj) ::   zqlw              ! long wave and sensible heat fluxes 
     757      REAL(wp), DIMENSION(jpi,jpj) ::   ztskk             ! skin temp. in Kelvin  
     758      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw              ! long wave and sensible heat fluxes       
    749759      REAL(wp), DIMENSION(jpi,jpj) ::   zqla              ! latent heat fluxes and evaporation 
    750       REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
    751760      !!--------------------------------------------------------------------- 
    752761      ! 
    753762      ! local scalars ( place there for vector optimisation purposes) 
    754       zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    755  
    756  
     763 
     764      ztskk(:,:) = ptsk(:,:) + rt0  ! => ptsk in Kelvin rather than Celsius 
     765       
    757766      ! ----------------------------------------------------------------------------- ! 
    758767      !     III    Net longwave radiative FLUX                                        ! 
     
    760769 
    761770      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 
    762       !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 
    763       zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux 
    764  
    765       !  Turbulent fluxes over ocean 
    766       ! ----------------------------- 
     771      !! (ztskk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 
     772      zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*ztskk(:,:)*ztskk(:,:)*ztskk(:,:)*ztskk(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux 
     773 
     774      !  Latent flux over ocean 
     775      ! ----------------------- 
    767776 
    768777      ! use scalar version of L_vap() for AGRIF compatibility 
    769778      DO jj = 1, jpj 
    770779         DO ji = 1, jpi 
    771             zqla(ji,jj) = -1._wp * L_vap( zst(ji,jj) ) * pevp(ji,jj)    ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 
     780            zqla(ji,jj) = - L_vap( ztskk(ji,jj) ) * pevp(ji,jj)    ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 
    772781         ENDDO 
    773782      ENDDO 
     
    788797      qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                   &   ! Downward Non Solar 
    789798         &     - psnow(:,:) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
    790          &     - pevp(:,:) * pst(:,:) * rcp                          &   ! remove evap heat content at SST !LB??? pst is Celsius !? 
     799         &     - pevp(:,:) * ptsk(:,:) * rcp                         &   ! remove evap heat content at SST 
    791800         &     + ( pprec(:,:) - psnow(:,:) ) * rn_pfac               &   ! add liquid precip heat content at Tair 
    792801         &     * ( ptair(:,:) - rt0 ) * rcp                          & 
     
    815824         CALL iom_put( "qsr_oce"  ,   qsr  )               ! output downward solar heat over the ocean 
    816825         CALL iom_put( "qt_oce"   ,   qns+qsr )            ! output total downward heat over the ocean 
    817       ENDIF 
    818       ! 
    819       IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    820          CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius 
    821          CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference... 
    822826      ENDIF 
    823827      ! 
     
    11051109 
    11061110      IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN 
    1107          ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) )  
     1111         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 
    11081112         IF( iom_use('evap_ao_cea'  ) )  CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average) 
    11091113         IF( iom_use('hflx_evap_cea') )  CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) )   ! heat flux from evap (cell average) 
     
    11141118      ENDIF 
    11151119      IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN 
    1116           WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) ;   ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 
    1117           ELSEWHERE                             ;   ztmp(:,:) = rcp * sst_m(:,:)     
    1118           ENDWHERE 
    1119           ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus )  
    1120           IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average) 
    1121           IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 
    1122           IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice) 
     1120         WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) 
     1121            ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 
     1122         ELSEWHERE 
     1123            ztmp(:,:) = rcp * sst_m(:,:) 
     1124         ENDWHERE 
     1125         ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus ) 
     1126         IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average) 
     1127         IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 
     1128         IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice) 
    11231129      ENDIF 
    11241130      ! 
Note: See TracChangeset for help on using the changeset viewer.