Ignore:
Timestamp:
2019-11-25T23:31:07+01:00 (10 months ago)
Author:
laurent
Message:

Syntax improvements and minor bug fixes…

File:
1 edited

Legend:

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

    r11841 r11962  
    142142      !!------------------------------------------------------------------- 
    143143      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 
    144          &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 
     144         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), tsk(jpi,jpj), STAT=sbc_blk_alloc ) 
    145145      ! 
    146146      CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) 
    147147      IF( sbc_blk_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' ) 
    148148   END FUNCTION sbc_blk_alloc 
    149  
    150    !LB: 
    151    INTEGER FUNCTION sbc_blk_cswl_alloc() 
    152       !!------------------------------------------------------------------- 
    153       !!             ***  ROUTINE sbc_blk_cswl_alloc *** 
    154       !!------------------------------------------------------------------- 
    155       !WRITE(numout,*) '*** LB: allocating tsk!' 
    156       ALLOCATE( tsk(jpi,jpj), STAT=sbc_blk_cswl_alloc ) 
    157       !WRITE(numout,*) '*** LB: done!' 
    158       CALL mpp_sum ( 'sbcblk', sbc_blk_cswl_alloc ) 
    159       IF( sbc_blk_cswl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_cswl_alloc: failed to allocate arrays' ) 
    160    END FUNCTION sbc_blk_cswl_alloc 
    161    !LB. 
    162149 
    163150 
     
    222209      !LB: 
    223210      !                             !** initialization of the cool-skin / warm-layer parametrization 
    224       IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    225          IF ( ln_NCAR ) CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' ) 
    226          !                       ! allocate array(s) for cool-skin/warm-layer param. 
    227          IF( sbc_blk_cswl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 
    228       END IF 
     211      IF( ln_NCAR .AND. (ln_skin_cs .OR. ln_skin_wl) ) & 
     212         & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' ) 
    229213      ! 
    230214      ioptio = 0 
     
    322306         ! 
    323307         WRITE(numout,*) 
    324          WRITE(numout,*) '      use cool-skin  parameterization (SSST)  ln_skin_cs  = ', ln_skin_cs !LB 
    325          WRITE(numout,*) '      use warm-layer parameterization (SSST)  ln_skin_wl  = ', ln_skin_wl !LB 
     308         WRITE(numout,*) '      use cool-skin  parameterization (SSST)  ln_skin_cs  = ', ln_skin_cs 
     309         WRITE(numout,*) '      use warm-layer parameterization (SSST)  ln_skin_wl  = ', ln_skin_wl 
    326310         ! 
    327          !LB: 
    328311         WRITE(numout,*) 
    329312         SELECT CASE( nhumi )              !* Print the choice of air humidity 
     
    332315         CASE( np_humi_rlh )   ;   WRITE(numout,*) '   ==>>>   air humidity is RELATIVE HUMIDITY     [%]' 
    333316         END SELECT 
    334          !LB. 
    335317         ! 
    336318      ENDIF 
     
    390372         ENDIF 
    391373         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1) 
    392          !LB: 
     374 
    393375         SELECT CASE( nhumi ) 
    394376         CASE( np_humi_sph ) 
     
    399381            qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !LB: 0.01 => RH is % percent in file 
    400382         END SELECT 
    401          !LB. 
     383 
    402384         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac 
    403385         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac 
     
    445427      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    446428      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
    447       REAL(wp), DIMENSION(jpi,jpj) ::   zqair             ! specific humidity     of air at z=rn_zqt [kg/kg] !LB 
     429      REAL(wp), DIMENSION(jpi,jpj) ::   zqair             ! specific humidity     of air at z=rn_zqt [kg/kg] 
    448430      !!--------------------------------------------------------------------- 
    449431      ! 
     
    510492         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity! 
    511493      CASE( np_humi_dpt ) 
    512          IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 
     494         !IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 
    513495         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    514496      CASE( np_humi_rlh ) 
    515          IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 
     497         !IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 
    516498         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file 
    517499      END SELECT 
     
    529511 
    530512         CASE( np_COARE_3p0 ) 
    531             IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p0" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm 
    532             CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.0 
    533                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    & 
     513            CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     514               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     & 
    534515               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
    535516 
    536517         CASE( np_COARE_3p6 ) 
    537             IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p6" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm 
    538             CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6 
    539                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    & 
     518            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     519               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     & 
    540520               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
    541521 
    542522         CASE( np_ECMWF     ) 
    543             IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_ecmwf" WITH CSWL options!!!, gdept_1d(1)=', gdept_1d(1) !LBrm 
    544             CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF 
    545                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   & 
     523            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
     524               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     & 
    546525               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
    547526 
     
    566545      ELSE !IF ( ln_skin_cs .OR. ln_skin_wl ) 
    567546 
    568  
    569547         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    570548            ! 
    571549         CASE( np_NCAR      ) 
    572             CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2 
     550            CALL turb_ncar    (      rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,                         & 
    573551               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    574552 
    575553         CASE( np_COARE_3p0 ) 
    576             IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p0" WITHOUT CSWL optional arrays!!!' !LBrm 
    577             CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.0 
     554            CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
    578555               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    579556 
    580557         CASE( np_COARE_3p6 ) 
    581             IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_coare3p6" WITHOUT CSWL optional arrays!!!' !LBrm 
    582             CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6 
     558            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
    583559               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    584560 
    585561         CASE( np_ECMWF     ) 
    586             IF (lwp) WRITE(numout,*) ' *** blk_oce => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!' !LBrm 
    587             CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF 
     562            CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
    588563               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    589564 
     
    596571      !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
    597572      !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef. 
    598        
     573 
    599574      IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 
    600575         !! If zu == zt, then ensuring once for all that: 
     
    602577         q_zu(:,:) = zqair(:,:) 
    603578      END IF 
    604        
     579 
    605580 
    606581      !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90 
    607582      ! ------------------------------------------------------------- 
    608        
    609       CALL BULK_FORMULA( 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        
     583 
     584      CALL BULK_FORMULA( rn_zu, zst(:,:), zsq(:,:), t_zu(:,:), q_zu(:,:), & 
     585         &               Cd_atm(:,:), Ch_atm(:,:), Ce_atm(:,:),           & 
     586         &               wndm(:,:), zU_zu(:,:), sf(jp_slp)%fnow(:,:,1),   & 
     587         &               taum(:,:), zqsb(:,:), zqla(:,:),                 & 
     588         &               pEvap=zevap(:,:), prhoa=rhoa(:,:) ) 
     589 
    614590      zqla(:,:)  =  zqla(:,:) * tmask(:,:,1) 
    615591      zqsb(:,:)  =  zqsb(:,:) * tmask(:,:,1) 
    616592      taum(:,:)  =  taum(:,:) * tmask(:,:,1) 
    617593      zevap(:,:) = zevap(:,:) * tmask(:,:,1) 
    618        
     594 
    619595      ! Tau i and j component on T-grid points, using array "Cd_atm" as a temporary array... 
    620596      Cd_atm = 0._wp 
    621597      WHERE ( wndm > 0._wp ) Cd_atm = taum / wndm 
    622598      zwnd_i = Cd_atm * zwnd_i 
    623       zwnd_j = Cd_atm * zwnd_j       
    624       !DO jj = 1, jpj             ! tau i and j component on T-grid points 
    625       !   DO ji = 1, jpi 
    626       !      zztmp = taum(ji,jj) / MAX( wndm(ji,jj) , 0.01_wp ) 
    627       !      zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
    628       !      zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
    629       !   END DO 
    630       !END DO 
     599      zwnd_j = Cd_atm * zwnd_j 
     600 
    631601      !                          ! add the HF tau contribution to the wind stress module 
    632602      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 
     
    652622      ! ----------------------------------------------------------------------------- ! 
    653623 
    654       !! 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.) 
     624      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 
     625      !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 
    655626      zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux 
    656627 
     
    848819         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity! 
    849820      CASE( np_humi_dpt ) 
    850          IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm 
     821         !IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm 
    851822         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    852823      CASE( np_humi_rlh ) 
    853          IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm 
     824         !IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm 
    854825         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file 
    855826      END SELECT 
Note: See TracChangeset for help on using the changeset viewer.