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 5407 for trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90 – NEMO

Ignore:
Timestamp:
2015-06-11T21:13:22+02:00 (9 years ago)
Author:
smasson
Message:

merge dev_r5218_CNRS17_coupling into the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r5385 r5407  
    4444   USE sbc_ice         ! Surface boundary condition: ice fields 
    4545   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 
    4653 
    4754   IMPLICIT NONE 
     
    4956 
    5057   PUBLIC   sbc_blk_core         ! routine called in sbcmod module 
    51    PUBLIC   blk_ice_core         ! 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 
    5262   PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module 
    5363 
     
    371381      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.) 
    372382         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
    373       qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar flux 
     383      ! 
     384      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar  
    374385         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip 
    375386         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
     
    379390         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) 
    380391      ! 
    381       CALL iom_put( "qlw_oce",   zqlw )                 ! output downward longwave heat over the ocean 
    382       CALL iom_put( "qsb_oce", - zqsb )                 ! output downward sensible heat over the ocean 
    383       CALL iom_put( "qla_oce", - zqla )                 ! output downward latent   heat over the ocean 
    384       CALL iom_put( "qhc_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean 
    385       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 
    386406      ! 
    387407      IF(ln_ctl) THEN 
     
    401421  
    402422    
    403    SUBROUTINE blk_ice_core(  pst   , pui   , pvi   , palb ,   & 
    404       &                      p_taui, p_tauj, p_qns , p_qsr,   & 
    405       &                      p_qla , p_dqns, p_dqla,          & 
    406       &                      p_tpr , p_spr ,                  & 
    407       &                      p_fr1 , p_fr2 , cd_grid, pdim  )  
    408       !!--------------------------------------------------------------------- 
    409       !!                     ***  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  *** 
    410427      !! 
    411428      !! ** Purpose :   provide the surface boundary condition over sea-ice 
    412429      !! 
    413       !! ** Method  :   compute momentum, heat and freshwater exchanged 
    414       !!                between atmosphere and sea-ice using CORE bulk 
    415       !!                formulea, ice variables and read atmmospheric fields. 
     430      !! ** Method  :   compute momentum using CORE bulk 
     431      !!                formulea, ice variables and read atmospheric fields. 
    416432      !!                NB: ice drag coefficient is assumed to be a constant 
    417       !!  
    418       !! caution : the net upward water flux has with mm/day unit 
    419       !!--------------------------------------------------------------------- 
    420       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin] 
    421       REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pui      ! ice surface velocity (i- and i- components      [m/s] 
    422       REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid) 
    423       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb     ! ice albedo (all skies)                            [%] 
    424       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_taui   ! i- & j-components of surface ice stress        [N/m2] 
    425       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid) 
    426       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2] 
    427       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2] 
    428       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2] 
    429       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2] 
    430       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2] 
    431       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s] 
    432       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s] 
    433       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%] 
    434       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%] 
    435       CHARACTER(len=1)          , INTENT(in   ) ::   cd_grid  ! ice grid ( C or B-grid) 
    436       INTEGER                   , INTENT(in   ) ::   pdim     ! number of ice categories 
    437       !! 
    438       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    439       INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays) 
    440       REAL(wp) ::   zst2, zst3 
    441       REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2, zcoef_dqlw, zcoef_dqla, zcoef_dqsb 
    442       REAL(wp) ::   zztmp                                        ! temporary variable 
    443       REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f                  ! relative wind module and components at F-point 
    444       REAL(wp) ::             zwndi_t , zwndj_t                  ! relative wind components at T-point 
    445       !! 
    446       REAL(wp), DIMENSION(:,:)  , POINTER ::   z_wnds_t          ! wind speed ( = | U10m - U_ice | ) at T-point 
    447       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice 
    448       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice 
    449       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice 
    450       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice 
    451       !!--------------------------------------------------------------------- 
    452       ! 
    453       IF( nn_timing == 1 )  CALL timing_start('blk_ice_core') 
    454       ! 
    455       CALL wrk_alloc( jpi,jpj, z_wnds_t ) 
    456       CALL wrk_alloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb )  
    457  
    458       ijpl  = pdim                            ! number of ice categories 
    459  
     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      ! 
    460442      ! local scalars ( place there for vector optimisation purposes) 
    461443      zcoef_wnorm  = rhoa * Cice 
    462444      zcoef_wnorm2 = rhoa * Cice * 0.5 
    463       zcoef_dqlw   = 4.0 * 0.95 * Stef 
    464       zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8) 
    465       zcoef_dqsb   = rhoa * cpa * Cice 
    466445 
    467446!!gm brutal.... 
    468       z_wnds_t(:,:) = 0.e0 
    469       p_taui  (:,:) = 0.e0 
    470       p_tauj  (:,:) = 0.e0 
     447      utau_ice  (:,:) = 0._wp 
     448      vtau_ice  (:,:) = 0._wp 
     449      wndm_ice  (:,:) = 0._wp 
    471450!!gm end 
    472451 
    473 #if defined key_lim3 
    474       tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)   ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init 
    475 #endif 
    476452      ! ----------------------------------------------------------------------------- ! 
    477453      !    Wind components and module relative to the moving ocean ( U10m - U_ice )   ! 
    478454      ! ----------------------------------------------------------------------------- ! 
    479       SELECT CASE( cd_grid ) 
     455      SELECT CASE( cp_ice_msh ) 
    480456      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
    481457         !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked) 
     
    484460               ! ... scalar wind at I-point (fld being at T-point) 
    485461               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
    486                   &              + 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) 
    487463               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    488                   &              + 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) 
    489465               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    490466               ! ... ice stress at I-point 
    491                p_taui(ji,jj) = zwnorm_f * zwndi_f 
    492                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 
    493469               ! ... scalar wind at T-point (fld being at T-point) 
    494                zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  pui(ji,jj+1) + pui(ji+1,jj+1)   & 
    495                   &                                                    + pui(ji,jj  ) + pui(ji+1,jj  )  ) 
    496                zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  pvi(ji,jj+1) + pvi(ji+1,jj+1)   & 
    497                   &                                                    + pvi(ji,jj  ) + pvi(ji+1,jj  )  ) 
    498                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) 
    499475            END DO 
    500476         END DO 
    501          CALL lbc_lnk( p_taui  , 'I', -1. ) 
    502          CALL lbc_lnk( p_tauj  , 'I', -1. ) 
    503          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. ) 
    504480         ! 
    505481      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    506482         DO jj = 2, jpj 
    507483            DO ji = fs_2, jpi   ! vect. opt. 
    508                zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pui(ji-1,jj  ) + pui(ji,jj) )  ) 
    509                zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pvi(ji  ,jj-1) + pvi(ji,jj) )  ) 
    510                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) 
    511487            END DO 
    512488         END DO 
    513489         DO jj = 2, jpjm1 
    514490            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    515                p_taui(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji+1,jj  ) + z_wnds_t(ji,jj) )                          & 
    516                   &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * pui(ji,jj) ) 
    517                p_tauj(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji,jj+1  ) + z_wnds_t(ji,jj) )                          & 
    518                   &          * ( 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) ) 
    519495            END DO 
    520496         END DO 
    521          CALL lbc_lnk( p_taui  , 'U', -1. ) 
    522          CALL lbc_lnk( p_tauj  , 'V', -1. ) 
    523          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. ) 
    524500         ! 
    525501      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 
    526548 
    527549      zztmp = 1. / ( 1. - albo ) 
    528550      !                                     ! ========================== ! 
    529       DO jl = 1, ijpl                       !  Loop over ice categories  ! 
     551      DO jl = 1, jpl                        !  Loop over ice categories  ! 
    530552         !                                  ! ========================== ! 
    531553         DO jj = 1 , jpj 
     
    534556               !      I   Radiative FLUXES   ! 
    535557               ! ----------------------------! 
    536                zst2 = pst(ji,jj,jl) * pst(ji,jj,jl) 
    537                zst3 = pst(ji,jj,jl) * zst2 
     558               zst2 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) 
     559               zst3 = ptsu(ji,jj,jl) * zst2 
    538560               ! Short Wave (sw) 
    539                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) 
    540562               ! Long  Wave (lw) 
    541                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) 
    542564               ! lw sensitivity 
    543565               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                                
     
    549571               ! ... turbulent heat fluxes 
    550572               ! Sensible Heat 
    551                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) ) 
    552574               ! Latent Heat 
    553                p_qla(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * z_wnds_t(ji,jj)   &                            
    554                   &                         * (  11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    555                ! Latent heat sensitivity for ice (Dqla/Dt) 
    556                IF( p_qla(ji,jj,jl) > 0._wp ) THEN 
    557                   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) ) 
    558580               ELSE 
    559                   p_dqla(ji,jj,jl) = 0._wp 
     581                  dqla_ice(ji,jj,jl) = 0._wp 
    560582               ENDIF 
    561583 
    562584               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    563                z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj) 
     585               z_dqsb(ji,jj,jl) = zcoef_dqsb * wndm_ice(ji,jj) 
    564586 
    565587               ! ----------------------------! 
     
    567589               ! ----------------------------! 
    568590               ! Downward Non Solar flux 
    569                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) 
    570592               ! Total non solar heat flux sensitivity for ice 
    571                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) ) 
    572594            END DO 
    573595            ! 
     
    576598      END DO 
    577599      ! 
     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      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing  
     616      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
     617      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
     618      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) 
     619 
     620      ! --- heat flux associated with emp --- ! 
     621      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap at sst 
     622         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair 
     623         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow) 
     624         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     625      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only) 
     626         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     627 
     628      ! --- total solar and non solar fluxes --- ! 
     629      qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:) 
     630      qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 
     631 
     632      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
     633      qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     634 
     635      CALL wrk_dealloc( jpi,jpj, zevap, zsnw )  
     636#endif 
     637 
    578638      !-------------------------------------------------------------------- 
    579639      ! FRACTIONs of net shortwave radiation which is not absorbed in the 
     
    581641      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    582642      ! 
    583       p_fr1(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    584       p_fr2(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
    585       ! 
    586       p_tpr(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s] 
    587       p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s] 
    588       CALL iom_put( 'snowpre', p_spr * 86400. )                  ! Snow precipitation 
    589       CALL iom_put( 'precip' , p_tpr * 86400. )                  ! Total precipitation 
     643      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     644      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     645      ! 
    590646      ! 
    591647      IF(ln_ctl) THEN 
    592          CALL prt_ctl(tab3d_1=p_qla   , clinfo1=' blk_ice_core: p_qla  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb  : ', kdim=ijpl) 
    593          CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw  : ', tab3d_2=p_dqla  , clinfo2=' p_dqla : ', kdim=ijpl) 
    594          CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw : ', kdim=ijpl) 
    595          CALL prt_ctl(tab3d_1=p_dqns  , clinfo1=' blk_ice_core: p_dqns : ', tab3d_2=p_qsr   , clinfo2=' p_qsr  : ', kdim=ijpl) 
    596          CALL prt_ctl(tab3d_1=pst     , clinfo1=' blk_ice_core: pst    : ', tab3d_2=p_qns   , clinfo2=' p_qns  : ', kdim=ijpl) 
    597          CALL prt_ctl(tab2d_1=p_tpr   , clinfo1=' blk_ice_core: p_tpr  : ', tab2d_2=p_spr   , clinfo2=' p_spr  : ') 
    598          CALL prt_ctl(tab2d_1=p_taui  , clinfo1=' blk_ice_core: p_taui : ', tab2d_2=p_tauj  , clinfo2=' p_tauj : ') 
    599          CALL prt_ctl(tab2d_1=z_wnds_t, clinfo1=' blk_ice_core: z_wnds_t : ') 
    600       ENDIF 
    601  
    602       CALL wrk_dealloc( jpi,jpj,   z_wnds_t ) 
    603       CALL wrk_dealloc( jpi,jpj,   pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    604       ! 
    605       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core') 
    606       ! 
    607    END SUBROUTINE blk_ice_core 
     648         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice_core: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl) 
     649         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl) 
     650         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl) 
     651         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice_core: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl) 
     652         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice_core: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl) 
     653         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice_core: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ') 
     654      ENDIF 
     655 
     656      CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
     657      ! 
     658      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx') 
     659       
     660   END SUBROUTINE blk_ice_core_flx 
     661#endif 
    608662 
    609663   SUBROUTINE turb_core_2z( zt, zu, sst, T_zt, q_sat, q_zt, dU,    & 
Note: See TracChangeset for help on using the changeset viewer.