Ignore:
Timestamp:
2015-07-15T17:46:12+02:00 (5 years ago)
Author:
andrewryan
Message:

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r5034 r5600  
    2222   !!   blk_oce_core    : computes momentum, heat and freshwater fluxes over ocean 
    2323   !!   blk_ice_core    : computes momentum, heat and freshwater fluxes over ice 
    24    !!   blk_bio_meanqsr : compute daily mean short wave radiation over the ocean 
    25    !!   blk_ice_meanqsr : compute daily mean short wave radiation over the ice 
    2624   !!   turb_core_2z    : Computes turbulent transfert coefficients 
    2725   !!   cd_neutral_10m  : Estimate of the neutral drag coefficient at 10m 
     
    4644   USE sbc_ice         ! Surface boundary condition: ice fields 
    4745   USE lib_fortran     ! to use key_nosignedzero 
     46#if defined key_lim3 
     47   USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b 
     48   USE limthd_dh       ! for CALL lim_thd_snwblow 
     49#elif defined key_lim2 
     50   USE ice_2, ONLY     : u_ice, v_ice 
     51   USE par_ice_2 
     52#endif 
    4853 
    4954   IMPLICIT NONE 
     
    5156 
    5257   PUBLIC   sbc_blk_core         ! routine called in sbcmod module 
    53    PUBLIC   blk_ice_core         ! routine called in sbc_ice_lim module 
    54    PUBLIC   blk_ice_meanqsr      ! routine called in sbc_ice_lim module 
     58#if defined key_lim2 || defined key_lim3 
     59   PUBLIC   blk_ice_core_tau     ! routine called in sbc_ice_lim module 
     60   PUBLIC   blk_ice_core_flx     ! routine called in sbc_ice_lim module 
     61#endif 
    5562   PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module 
    5663 
     
    195202      !                                            ! compute the surface ocean fluxes using CORE bulk formulea 
    196203      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce_core( kt, sf, sst_m, ssu_m, ssv_m ) 
    197  
    198       ! If diurnal cycle is activated, compute a daily mean short waves flux for biogeochemistery 
    199       IF( ltrcdm2dc )   CALL blk_bio_meanqsr 
    200204 
    201205#if defined key_cice 
     
    302306      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
    303307      ENDIF 
     308 
    304309      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    305310      ! ----------------------------------------------------------------------------- ! 
     
    376381      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.) 
    377382         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    378       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar flux 
     383      ! 
     384      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar  
    379385         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip 
    380386         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
     
    384390         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) 
    385391      ! 
    386       CALL iom_put( "qlw_oce",   zqlw )                 ! output downward longwave heat over the ocean 
    387       CALL iom_put( "qsb_oce", - zqsb )                 ! output downward sensible heat over the ocean 
    388       CALL iom_put( "qla_oce", - zqla )                 ! output downward latent   heat over the ocean 
    389       CALL iom_put( "qhc_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean 
    390       CALL iom_put( "qns_oce",   qns  )                 ! output downward non solar heat over the ocean 
     392#if defined key_lim3 
     393      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by LIM3) 
     394      qsr_oce(:,:) = qsr(:,:) 
     395#endif 
     396      ! 
     397      IF ( nn_ice == 0 ) THEN 
     398         CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean 
     399         CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean 
     400         CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean 
     401         CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean 
     402         CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean 
     403         CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean 
     404         CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean 
     405         tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac   ! output total precipitation [kg/m2/s] 
     406         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac   ! output solid precipitation [kg/m2/s] 
     407         CALL iom_put( 'snowpre', sprecip * 86400. )        ! Snow 
     408         CALL iom_put( 'precip' , tprecip * 86400. )        ! Total precipitation 
     409      ENDIF 
    391410      ! 
    392411      IF(ln_ctl) THEN 
     
    406425  
    407426    
    408    SUBROUTINE blk_ice_core(  pst   , pui   , pvi   , palb ,   & 
    409       &                      p_taui, p_tauj, p_qns , p_qsr,   & 
    410       &                      p_qla , p_dqns, p_dqla,          & 
    411       &                      p_tpr , p_spr ,                  & 
    412       &                      p_fr1 , p_fr2 , cd_grid, pdim  )  
    413       !!--------------------------------------------------------------------- 
    414       !!                     ***  ROUTINE blk_ice_core  *** 
     427#if defined key_lim2 || defined key_lim3 
     428   SUBROUTINE blk_ice_core_tau 
     429      !!--------------------------------------------------------------------- 
     430      !!                     ***  ROUTINE blk_ice_core_tau  *** 
    415431      !! 
    416432      !! ** Purpose :   provide the surface boundary condition over sea-ice 
    417433      !! 
    418       !! ** Method  :   compute momentum, heat and freshwater exchanged 
    419       !!                between atmosphere and sea-ice using CORE bulk 
    420       !!                formulea, ice variables and read atmmospheric fields. 
     434      !! ** Method  :   compute momentum using CORE bulk 
     435      !!                formulea, ice variables and read atmospheric fields. 
    421436      !!                NB: ice drag coefficient is assumed to be a constant 
    422       !!  
    423       !! caution : the net upward water flux has with mm/day unit 
    424       !!--------------------------------------------------------------------- 
    425       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin] 
    426       REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pui      ! ice surface velocity (i- and i- components      [m/s] 
    427       REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid) 
    428       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb     ! ice albedo (all skies)                            [%] 
    429       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_taui   ! i- & j-components of surface ice stress        [N/m2] 
    430       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid) 
    431       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2] 
    432       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2] 
    433       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2] 
    434       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2] 
    435       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2] 
    436       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s] 
    437       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s] 
    438       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%] 
    439       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%] 
    440       CHARACTER(len=1)          , INTENT(in   ) ::   cd_grid  ! ice grid ( C or B-grid) 
    441       INTEGER                   , INTENT(in   ) ::   pdim     ! number of ice categories 
    442       !! 
    443       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    444       INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays) 
    445       REAL(wp) ::   zst2, zst3 
    446       REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2, zcoef_dqlw, zcoef_dqla, zcoef_dqsb 
    447       REAL(wp) ::   zztmp                                        ! temporary variable 
    448       REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f                  ! relative wind module and components at F-point 
    449       REAL(wp) ::             zwndi_t , zwndj_t                  ! relative wind components at T-point 
    450       !! 
    451       REAL(wp), DIMENSION(:,:)  , POINTER ::   z_wnds_t          ! wind speed ( = | U10m - U_ice | ) at T-point 
    452       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice 
    453       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice 
    454       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice 
    455       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice 
    456       !!--------------------------------------------------------------------- 
    457       ! 
    458       IF( nn_timing == 1 )  CALL timing_start('blk_ice_core') 
    459       ! 
    460       CALL wrk_alloc( jpi,jpj, z_wnds_t ) 
    461       CALL wrk_alloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb )  
    462  
    463       ijpl  = pdim                            ! number of ice categories 
    464  
     437      !!--------------------------------------------------------------------- 
     438      INTEGER  ::   ji, jj    ! dummy loop indices 
     439      REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2 
     440      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point 
     441      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point 
     442      !!--------------------------------------------------------------------- 
     443      ! 
     444      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_tau') 
     445      ! 
    465446      ! local scalars ( place there for vector optimisation purposes) 
    466447      zcoef_wnorm  = rhoa * Cice 
    467448      zcoef_wnorm2 = rhoa * Cice * 0.5 
    468       zcoef_dqlw   = 4.0 * 0.95 * Stef 
    469       zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
    470       zcoef_dqsb   = rhoa * cpa * Cice 
    471449 
    472450!!gm brutal.... 
    473       z_wnds_t(:,:) = 0.e0 
    474       p_taui  (:,:) = 0.e0 
    475       p_tauj  (:,:) = 0.e0 
     451      utau_ice  (:,:) = 0._wp 
     452      vtau_ice  (:,:) = 0._wp 
     453      wndm_ice  (:,:) = 0._wp 
    476454!!gm end 
    477455 
    478 #if defined key_lim3 
    479       tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)   ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init 
    480 #endif 
    481456      ! ----------------------------------------------------------------------------- ! 
    482457      !    Wind components and module relative to the moving ocean ( U10m - U_ice )   ! 
    483458      ! ----------------------------------------------------------------------------- ! 
    484       SELECT CASE( cd_grid ) 
     459      SELECT CASE( cp_ice_msh ) 
    485460      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
    486461         !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked) 
     
    489464               ! ... scalar wind at I-point (fld being at T-point) 
    490465               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
    491                   &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * pui(ji,jj) 
     466                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
    492467               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    493                   &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * pvi(ji,jj) 
     468                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    494469               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    495470               ! ... ice stress at I-point 
    496                p_taui(ji,jj) = zwnorm_f * zwndi_f 
    497                p_tauj(ji,jj) = zwnorm_f * zwndj_f 
     471               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
     472               vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
    498473               ! ... scalar wind at T-point (fld being at T-point) 
    499                zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  pui(ji,jj+1) + pui(ji+1,jj+1)   & 
    500                   &                                                    + pui(ji,jj  ) + pui(ji+1,jj  )  ) 
    501                zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  pvi(ji,jj+1) + pvi(ji+1,jj+1)   & 
    502                   &                                                    + pvi(ji,jj  ) + pvi(ji+1,jj  )  ) 
    503                z_wnds_t(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     474               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
     475                  &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
     476               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
     477                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
     478               wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    504479            END DO 
    505480         END DO 
    506          CALL lbc_lnk( p_taui  , 'I', -1. ) 
    507          CALL lbc_lnk( p_tauj  , 'I', -1. ) 
    508          CALL lbc_lnk( z_wnds_t, 'T',  1. ) 
     481         CALL lbc_lnk( utau_ice, 'I', -1. ) 
     482         CALL lbc_lnk( vtau_ice, 'I', -1. ) 
     483         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    509484         ! 
    510485      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    511486         DO jj = 2, jpj 
    512487            DO ji = fs_2, jpi   ! vect. opt. 
    513                zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pui(ji-1,jj  ) + pui(ji,jj) )  ) 
    514                zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pvi(ji  ,jj-1) + pvi(ji,jj) )  ) 
    515                z_wnds_t(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     488               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
     489               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
     490               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    516491            END DO 
    517492         END DO 
    518493         DO jj = 2, jpjm1 
    519494            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    520                p_taui(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji+1,jj  ) + z_wnds_t(ji,jj) )                          & 
    521                   &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * pui(ji,jj) ) 
    522                p_tauj(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji,jj+1  ) + z_wnds_t(ji,jj) )                          & 
    523                   &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * pvi(ji,jj) ) 
     495               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     496                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
     497               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     498                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    524499            END DO 
    525500         END DO 
    526          CALL lbc_lnk( p_taui  , 'U', -1. ) 
    527          CALL lbc_lnk( p_tauj  , 'V', -1. ) 
    528          CALL lbc_lnk( z_wnds_t, 'T',  1. ) 
     501         CALL lbc_lnk( utau_ice, 'U', -1. ) 
     502         CALL lbc_lnk( vtau_ice, 'V', -1. ) 
     503         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    529504         ! 
    530505      END SELECT 
     506 
     507      IF(ln_ctl) THEN 
     508         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice_core: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ') 
     509         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice_core: wndm_ice : ') 
     510      ENDIF 
     511 
     512      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau') 
     513       
     514   END SUBROUTINE blk_ice_core_tau 
     515 
     516 
     517   SUBROUTINE blk_ice_core_flx( ptsu, palb ) 
     518      !!--------------------------------------------------------------------- 
     519      !!                     ***  ROUTINE blk_ice_core_flx  *** 
     520      !! 
     521      !! ** Purpose :   provide the surface boundary condition over sea-ice 
     522      !! 
     523      !! ** Method  :   compute heat and freshwater exchanged 
     524      !!                between atmosphere and sea-ice using CORE bulk 
     525      !!                formulea, ice variables and read atmmospheric fields. 
     526      !!  
     527      !! caution : the net upward water flux has with mm/day unit 
     528      !!--------------------------------------------------------------------- 
     529      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu          ! sea ice surface temperature 
     530      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb          ! ice albedo (all skies) 
     531      !! 
     532      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     533      REAL(wp) ::   zst2, zst3 
     534      REAL(wp) ::   zcoef_dqlw, zcoef_dqla, zcoef_dqsb 
     535      REAL(wp) ::   zztmp, z1_lsub                               ! temporary variable 
     536      !! 
     537      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice 
     538      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice 
     539      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice 
     540      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice 
     541      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw       ! evaporation and snw distribution after wind blowing (LIM3) 
     542      !!--------------------------------------------------------------------- 
     543      ! 
     544      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx') 
     545      ! 
     546      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )  
     547 
     548      ! local scalars ( place there for vector optimisation purposes) 
     549      zcoef_dqlw   = 4.0 * 0.95 * Stef 
     550      zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
     551      zcoef_dqsb   = rhoa * cpa * Cice 
    531552 
    532553      zztmp = 1. / ( 1. - albo ) 
    533554      !                                     ! ========================== ! 
    534       DO jl = 1, ijpl                       !  Loop over ice categories  ! 
     555      DO jl = 1, jpl                        !  Loop over ice categories  ! 
    535556         !                                  ! ========================== ! 
    536557         DO jj = 1 , jpj 
     
    539560               !      I   Radiative FLUXES   ! 
    540561               ! ----------------------------! 
    541                zst2 = pst(ji,jj,jl) * pst(ji,jj,jl) 
    542                zst3 = pst(ji,jj,jl) * zst2 
     562               zst2 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) 
     563               zst3 = ptsu(ji,jj,jl) * zst2 
    543564               ! Short Wave (sw) 
    544                p_qsr(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
     565               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
    545566               ! Long  Wave (lw) 
    546                z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
     567               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    547568               ! lw sensitivity 
    548569               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                                
     
    554575               ! ... turbulent heat fluxes 
    555576               ! Sensible Heat 
    556                z_qsb(ji,jj,jl) = rhoa * cpa * Cice * z_wnds_t(ji,jj) * ( pst(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
     577               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    557578               ! Latent Heat 
    558                p_qla(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * z_wnds_t(ji,jj)   &                            
    559                   &                         * (  11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    560                ! Latent heat sensitivity for ice (Dqla/Dt) 
    561                IF( p_qla(ji,jj,jl) > 0._wp ) THEN 
    562                   p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) ) 
     579               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                            
     580                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
     581              ! Latent heat sensitivity for ice (Dqla/Dt) 
     582               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
     583                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 
    563584               ELSE 
    564                   p_dqla(ji,jj,jl) = 0._wp 
     585                  dqla_ice(ji,jj,jl) = 0._wp 
    565586               ENDIF 
    566587 
    567588               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    568                z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj) 
     589               z_dqsb(ji,jj,jl) = zcoef_dqsb * wndm_ice(ji,jj) 
    569590 
    570591               ! ----------------------------! 
     
    572593               ! ----------------------------! 
    573594               ! Downward Non Solar flux 
    574                p_qns (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - p_qla (ji,jj,jl) 
     595               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) 
    575596               ! Total non solar heat flux sensitivity for ice 
    576                p_dqns(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + p_dqla(ji,jj,jl) ) 
     597               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) 
    577598            END DO 
    578599            ! 
     
    581602      END DO 
    582603      ! 
     604      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s] 
     605      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s] 
     606      CALL iom_put( 'snowpre', sprecip * 86400. )                  ! Snow precipitation 
     607      CALL iom_put( 'precip' , tprecip * 86400. )                  ! Total precipitation 
     608 
     609#if defined  key_lim3 
     610      CALL wrk_alloc( jpi,jpj, zevap, zsnw )  
     611 
     612      ! --- evaporation --- ! 
     613      z1_lsub = 1._wp / Lsub 
     614      evap_ice (:,:,:) = qla_ice (:,:,:) * z1_lsub ! sublimation 
     615      devap_ice(:,:,:) = dqla_ice(:,:,:) * z1_lsub 
     616      zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean 
     617 
     618      ! --- evaporation minus precipitation --- ! 
     619      zsnw(:,:) = 0._wp 
     620      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing  
     621      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
     622      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
     623      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) 
     624 
     625      ! --- heat flux associated with emp --- ! 
     626      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap at sst 
     627         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair 
     628         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow) 
     629         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     630      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only) 
     631         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     632 
     633      ! --- total solar and non solar fluxes --- ! 
     634      qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:) 
     635      qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 
     636 
     637      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
     638      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     639 
     640      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  
     641#endif 
     642 
    583643      !-------------------------------------------------------------------- 
    584644      ! FRACTIONs of net shortwave radiation which is not absorbed in the 
     
    586646      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    587647      ! 
    588       p_fr1(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    589       p_fr2(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
    590       ! 
    591       p_tpr(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s] 
    592       p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s] 
    593       CALL iom_put( 'snowpre', p_spr * 86400. )                  ! Snow precipitation 
    594       CALL iom_put( 'precip' , p_tpr * 86400. )                  ! Total precipitation 
     648      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     649      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     650      ! 
    595651      ! 
    596652      IF(ln_ctl) THEN 
    597          CALL prt_ctl(tab3d_1=p_qla   , clinfo1=' blk_ice_core: p_qla  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb  : ', kdim=ijpl) 
    598          CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw  : ', tab3d_2=p_dqla  , clinfo2=' p_dqla : ', kdim=ijpl) 
    599          CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw : ', kdim=ijpl) 
    600          CALL prt_ctl(tab3d_1=p_dqns  , clinfo1=' blk_ice_core: p_dqns : ', tab3d_2=p_qsr   , clinfo2=' p_qsr  : ', kdim=ijpl) 
    601          CALL prt_ctl(tab3d_1=pst     , clinfo1=' blk_ice_core: pst    : ', tab3d_2=p_qns   , clinfo2=' p_qns  : ', kdim=ijpl) 
    602          CALL prt_ctl(tab2d_1=p_tpr   , clinfo1=' blk_ice_core: p_tpr  : ', tab2d_2=p_spr   , clinfo2=' p_spr  : ') 
    603          CALL prt_ctl(tab2d_1=p_taui  , clinfo1=' blk_ice_core: p_taui : ', tab2d_2=p_tauj  , clinfo2=' p_tauj : ') 
    604          CALL prt_ctl(tab2d_1=z_wnds_t, clinfo1=' blk_ice_core: z_wnds_t : ') 
    605       ENDIF 
    606  
    607       CALL wrk_dealloc( jpi,jpj,   z_wnds_t ) 
    608       CALL wrk_dealloc( jpi,jpj,   pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    609       ! 
    610       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core') 
    611       ! 
    612    END SUBROUTINE blk_ice_core 
    613  
    614  
    615    SUBROUTINE blk_bio_meanqsr 
    616       !!--------------------------------------------------------------------- 
    617       !!                     ***  ROUTINE blk_bio_meanqsr 
    618       !!                      
    619       !! ** Purpose :   provide daily qsr_mean for PISCES when 
    620       !!                analytic diurnal cycle is applied in physic 
    621       !!                 
    622       !! ** Method  :   add part where there is no ice 
    623       !!  
    624       !!--------------------------------------------------------------------- 
    625       IF( nn_timing == 1 )  CALL timing_start('blk_bio_meanqsr') 
    626       ! 
    627       qsr_mean(:,:) = (1. - albo ) *  sf(jp_qsr)%fnow(:,:,1) 
    628       ! 
    629       IF( nn_timing == 1 )  CALL timing_stop('blk_bio_meanqsr') 
    630       ! 
    631    END SUBROUTINE blk_bio_meanqsr 
    632   
    633   
    634    SUBROUTINE blk_ice_meanqsr( palb, p_qsr_mean, pdim ) 
    635       !!--------------------------------------------------------------------- 
    636       !! 
    637       !! ** Purpose :   provide the daily qsr_mean over sea_ice for PISCES when 
    638       !!                analytic diurnal cycle is applied in physic 
    639       !! 
    640       !! ** Method  :   compute qsr 
    641       !!  
    642       !!--------------------------------------------------------------------- 
    643       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb       ! ice albedo (clear sky) (alb_ice_cs)               [%] 
    644       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qsr_mean !     solar heat flux over ice (T-point)         [W/m2] 
    645       INTEGER                   , INTENT(in   ) ::   pdim       ! number of ice categories 
    646       ! 
    647       INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays) 
    648       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    649       REAL(wp) ::   zztmp         ! temporary variable 
    650       !!--------------------------------------------------------------------- 
    651       IF( nn_timing == 1 )  CALL timing_start('blk_ice_meanqsr') 
    652       ! 
    653       ijpl  = pdim                            ! number of ice categories 
    654       zztmp = 1. / ( 1. - albo ) 
    655       !                                     ! ========================== ! 
    656       DO jl = 1, ijpl                       !  Loop over ice categories  ! 
    657          !                                  ! ========================== ! 
    658          DO jj = 1 , jpj 
    659             DO ji = 1, jpi 
    660                   p_qsr_mean(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr_mean(ji,jj) 
    661             END DO 
    662          END DO 
    663       END DO 
    664       ! 
    665       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_meanqsr') 
    666       ! 
    667    END SUBROUTINE blk_ice_meanqsr   
    668  
     653         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice_core: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl) 
     654         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl) 
     655         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl) 
     656         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice_core: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl) 
     657         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice_core: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl) 
     658         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice_core: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ') 
     659      ENDIF 
     660 
     661      CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
     662      ! 
     663      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx') 
     664       
     665   END SUBROUTINE blk_ice_core_flx 
     666#endif 
    669667 
    670668   SUBROUTINE turb_core_2z( zt, zu, sst, T_zt, q_sat, q_zt, dU,    & 
     
    848846      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1   
    849847      cd_neutral_10m = 1.e-3 * ( & 
    850          &       (rgt33 + 1._wp)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. 
     848         &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. 
    851849         &      + rgt33         *      2.34   )                                                    ! zw10 >= 33. 
    852850      ! 
Note: See TracChangeset for help on using the changeset viewer.