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 5709 for branches/UKMO/DEV_r5107_dynvor_updates/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 – NEMO

Ignore:
Timestamp:
2015-08-26T16:14:53+02:00 (9 years ago)
Author:
davestorkey
Message:

Updating UKMO/DEV_r5107_dynvor_updates branch to be relative to revision 5518 of the trunk
(= branching point for NEMO 3.6_STABLE).

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/DEV_r5107_dynvor_updates/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r5256 r5709  
    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      ENDIF 
    391406      ! 
    392407      IF(ln_ctl) THEN 
     
    406421  
    407422    
    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  *** 
     423#if defined key_lim2 || defined key_lim3 
     424   SUBROUTINE blk_ice_core_tau 
     425      !!--------------------------------------------------------------------- 
     426      !!                     ***  ROUTINE blk_ice_core_tau  *** 
    415427      !! 
    416428      !! ** Purpose :   provide the surface boundary condition over sea-ice 
    417429      !! 
    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. 
     430      !! ** Method  :   compute momentum using CORE bulk 
     431      !!                formulea, ice variables and read atmospheric fields. 
    421432      !!                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  
     433      !!--------------------------------------------------------------------- 
     434      INTEGER  ::   ji, jj    ! dummy loop indices 
     435      REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2 
     436      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point 
     437      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point 
     438      !!--------------------------------------------------------------------- 
     439      ! 
     440      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_tau') 
     441      ! 
    465442      ! local scalars ( place there for vector optimisation purposes) 
    466443      zcoef_wnorm  = rhoa * Cice 
    467444      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 
    471445 
    472446!!gm brutal.... 
    473       z_wnds_t(:,:) = 0.e0 
    474       p_taui  (:,:) = 0.e0 
    475       p_tauj  (:,:) = 0.e0 
     447      utau_ice  (:,:) = 0._wp 
     448      vtau_ice  (:,:) = 0._wp 
     449      wndm_ice  (:,:) = 0._wp 
    476450!!gm end 
    477451 
    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 
    481452      ! ----------------------------------------------------------------------------- ! 
    482453      !    Wind components and module relative to the moving ocean ( U10m - U_ice )   ! 
    483454      ! ----------------------------------------------------------------------------- ! 
    484       SELECT CASE( cd_grid ) 
     455      SELECT CASE( cp_ice_msh ) 
    485456      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
    486457         !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked) 
     
    489460               ! ... scalar wind at I-point (fld being at T-point) 
    490461               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) 
     462                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj) 
    492463               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) 
     464                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    494465               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    495466               ! ... ice stress at I-point 
    496                p_taui(ji,jj) = zwnorm_f * zwndi_f 
    497                p_tauj(ji,jj) = zwnorm_f * zwndj_f 
     467               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
     468               vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
    498469               ! ... 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) 
     470               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
     471                  &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
     472               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
     473                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
     474               wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    504475            END DO 
    505476         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. ) 
     477         CALL lbc_lnk( utau_ice, 'I', -1. ) 
     478         CALL lbc_lnk( vtau_ice, 'I', -1. ) 
     479         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    509480         ! 
    510481      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    511482         DO jj = 2, jpj 
    512483            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) 
     484               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
     485               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
     486               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    516487            END DO 
    517488         END DO 
    518489         DO jj = 2, jpjm1 
    519490            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) ) 
     491               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     492                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
     493               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     494                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    524495            END DO 
    525496         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. ) 
     497         CALL lbc_lnk( utau_ice, 'U', -1. ) 
     498         CALL lbc_lnk( vtau_ice, 'V', -1. ) 
     499         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    529500         ! 
    530501      END SELECT 
     502 
     503      IF(ln_ctl) THEN 
     504         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice_core: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ') 
     505         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice_core: wndm_ice : ') 
     506      ENDIF 
     507 
     508      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau') 
     509       
     510   END SUBROUTINE blk_ice_core_tau 
     511 
     512 
     513   SUBROUTINE blk_ice_core_flx( ptsu, palb ) 
     514      !!--------------------------------------------------------------------- 
     515      !!                     ***  ROUTINE blk_ice_core_flx  *** 
     516      !! 
     517      !! ** Purpose :   provide the surface boundary condition over sea-ice 
     518      !! 
     519      !! ** Method  :   compute heat and freshwater exchanged 
     520      !!                between atmosphere and sea-ice using CORE bulk 
     521      !!                formulea, ice variables and read atmmospheric fields. 
     522      !!  
     523      !! caution : the net upward water flux has with mm/day unit 
     524      !!--------------------------------------------------------------------- 
     525      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu          ! sea ice surface temperature 
     526      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb          ! ice albedo (all skies) 
     527      !! 
     528      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     529      REAL(wp) ::   zst2, zst3 
     530      REAL(wp) ::   zcoef_dqlw, zcoef_dqla, zcoef_dqsb 
     531      REAL(wp) ::   zztmp, z1_lsub                               ! temporary variable 
     532      !! 
     533      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice 
     534      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice 
     535      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice 
     536      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice 
     537      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw       ! evaporation and snw distribution after wind blowing (LIM3) 
     538      !!--------------------------------------------------------------------- 
     539      ! 
     540      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx') 
     541      ! 
     542      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )  
     543 
     544      ! local scalars ( place there for vector optimisation purposes) 
     545      zcoef_dqlw   = 4.0 * 0.95 * Stef 
     546      zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
     547      zcoef_dqsb   = rhoa * cpa * Cice 
    531548 
    532549      zztmp = 1. / ( 1. - albo ) 
    533550      !                                     ! ========================== ! 
    534       DO jl = 1, ijpl                       !  Loop over ice categories  ! 
     551      DO jl = 1, jpl                        !  Loop over ice categories  ! 
    535552         !                                  ! ========================== ! 
    536553         DO jj = 1 , jpj 
     
    539556               !      I   Radiative FLUXES   ! 
    540557               ! ----------------------------! 
    541                zst2 = pst(ji,jj,jl) * pst(ji,jj,jl) 
    542                zst3 = pst(ji,jj,jl) * zst2 
     558               zst2 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) 
     559               zst3 = ptsu(ji,jj,jl) * zst2 
    543560               ! Short Wave (sw) 
    544                p_qsr(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
     561               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
    545562               ! 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) 
     563               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    547564               ! lw sensitivity 
    548565               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                                
     
    554571               ! ... turbulent heat fluxes 
    555572               ! 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) ) 
     573               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    557574               ! 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) ) 
     575               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                            
     576                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
     577              ! Latent heat sensitivity for ice (Dqla/Dt) 
     578               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
     579                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 
    563580               ELSE 
    564                   p_dqla(ji,jj,jl) = 0._wp 
     581                  dqla_ice(ji,jj,jl) = 0._wp 
    565582               ENDIF 
    566583 
    567584               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    568                z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj) 
     585               z_dqsb(ji,jj,jl) = zcoef_dqsb * wndm_ice(ji,jj) 
    569586 
    570587               ! ----------------------------! 
     
    572589               ! ----------------------------! 
    573590               ! 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) 
     591               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) 
    575592               ! 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) ) 
     593               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) 
    577594            END DO 
    578595            ! 
     
    581598      END DO 
    582599      ! 
     600      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s] 
     601      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s] 
     602      CALL iom_put( 'snowpre', sprecip * 86400. )                  ! Snow precipitation 
     603      CALL iom_put( 'precip' , tprecip * 86400. )                  ! Total precipitation 
     604 
     605#if defined  key_lim3 
     606      CALL wrk_alloc( jpi,jpj, zevap, zsnw )  
     607 
     608      ! --- evaporation --- ! 
     609      z1_lsub = 1._wp / Lsub 
     610      evap_ice (:,:,:) = qla_ice (:,:,:) * z1_lsub ! sublimation 
     611      devap_ice(:,:,:) = dqla_ice(:,:,:) * z1_lsub 
     612      zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean 
     613 
     614      ! --- evaporation minus precipitation --- ! 
     615      zsnw(:,:) = 0._wp 
     616      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing  
     617      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
     618      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
     619      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) 
     620 
     621      ! --- heat flux associated with emp --- ! 
     622      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap at sst 
     623         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair 
     624         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow) 
     625         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     626      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only) 
     627         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     628 
     629      ! --- total solar and non solar fluxes --- ! 
     630      qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:) 
     631      qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 
     632 
     633      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
     634      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     635 
     636      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  
     637#endif 
     638 
    583639      !-------------------------------------------------------------------- 
    584640      ! FRACTIONs of net shortwave radiation which is not absorbed in the 
     
    586642      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    587643      ! 
    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 
     644      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     645      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     646      ! 
    595647      ! 
    596648      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  
     649         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice_core: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl) 
     650         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl) 
     651         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl) 
     652         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice_core: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl) 
     653         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice_core: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl) 
     654         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice_core: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ') 
     655      ENDIF 
     656 
     657      CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
     658      ! 
     659      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx') 
     660       
     661   END SUBROUTINE blk_ice_core_flx 
     662#endif 
    669663 
    670664   SUBROUTINE turb_core_2z( zt, zu, sst, T_zt, q_sat, q_zt, dU,    & 
Note: See TracChangeset for help on using the changeset viewer.