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 5038 for branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 – NEMO

Ignore:
Timestamp:
2015-01-20T15:26:13+01:00 (9 years ago)
Author:
jamesharle
Message:

Merging branch with HEAD of the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r4792 r5038  
    1212   !!            3.4  ! 2011-01  (A Porter)  dynamical allocation 
    1313   !!             -   ! 2012-10  (C. Rousset)  add lim_diahsb 
     14   !!            3.6  ! 2014-07  (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface 
    1415   !!---------------------------------------------------------------------- 
    1516#if defined key_lim3 
     
    8081   !!---------------------------------------------------------------------- 
    8182CONTAINS 
    82  
    83    FUNCTION fice_cell_ave ( ptab) 
    84       !!-------------------------------------------------------------------------- 
    85       !! * Compute average over categories, for grid cell (ice covered and free ocean) 
    86       !!-------------------------------------------------------------------------- 
    87       REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave 
    88       REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab 
    89       INTEGER :: jl ! Dummy loop index 
    90        
    91       fice_cell_ave (:,:) = 0.0_wp 
    92        
    93       DO jl = 1, jpl 
    94          fice_cell_ave (:,:) = fice_cell_ave (:,:) & 
    95             &                  + a_i (:,:,jl) * ptab (:,:,jl) 
    96       END DO 
    97        
    98    END FUNCTION fice_cell_ave 
    99     
    100    FUNCTION fice_ice_ave ( ptab) 
    101       !!-------------------------------------------------------------------------- 
    102       !! * Compute average over categories, for ice covered part of grid cell 
    103       !!-------------------------------------------------------------------------- 
    104       REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave 
    105       REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab 
    106  
    107       fice_ice_ave (:,:) = 0.0_wp 
    108       WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:) 
    109  
    110    END FUNCTION fice_ice_ave 
    11183 
    11284   !!====================================================================== 
     
    133105      !!--------------------------------------------------------------------- 
    134106      INTEGER, INTENT(in) ::   kt      ! ocean time step 
    135       INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE) 
     107      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE, =5 COUPLED) 
    136108      !! 
    137109      INTEGER  ::   ji, jj, jl, jk      ! dummy loop index 
    138110      REAL(wp) ::   zcoef   ! local scalar 
    139       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky 
    140       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice (for coupled) 
    141  
    142       REAL(wp), POINTER, DIMENSION(:,:) :: zalb_ice_all    ! Mean albedo over all categories 
    143       REAL(wp), POINTER, DIMENSION(:,:) :: ztem_ice_all    ! Mean temperature over all categories 
    144        
    145       REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_ice_all   ! Mean solar heat flux over all categories 
    146       REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_ice_all   ! Mean non solar heat flux over all categories 
    147       REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_ice_all   ! Mean latent heat flux over all categories 
    148       REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories 
    149       REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all  ! Mean d(qla)/dT over all categories 
    150       REAL(wp) ::   ztmelts           ! clem 2014: for HC diags 
    151       REAL(wp) ::   epsi20 = 1.e-20   ! 
     111      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky 
     112      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled) 
    152113      !!---------------------------------------------------------------------- 
    153114 
    154       !- O.M. : why do we allocate all these arrays even when MOD( kt-1, nn_fsbc ) /= 0 ????? 
    155  
    156115      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
    157  
    158       CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice ) 
    159  
    160       IF( lk_cpl ) THEN 
    161          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
    162             &   CALL wrk_alloc( jpi, jpj, ztem_ice_all , zalb_ice_all  , z_qsr_ice_all, z_qns_ice_all,   & 
    163             &                             z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
    164       ENDIF 
    165116 
    166117      IF( kt == nit000 ) THEN 
     
    183134         !                                           !----------------! 
    184135         ! 
    185          u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point 
    186          v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean) 
    187  
    188          ! masked sea surface freezing temperature [Kelvin] 
    189          t_bo(:,:) = ( tfreez( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) ) 
    190  
    191          CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo 
    192  
     136         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)                     ! mean surface ocean current at ice velocity point 
     137         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)                    ! (C-grid dynamics :  U- & V-points as the ocean) 
     138         ! 
     139         t_bo(:,:) = ( eos_fzp( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) )  ! masked sea surface freezing temperature [Kelvin] 
     140         !                                                                                  ! (set to rt0 over land) 
     141         !                                           ! Ice albedo 
     142         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )       
     143 
     144         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
     145 
     146         SELECT CASE( kblk ) 
     147         CASE( jp_core , jp_cpl )   ! CORE and COUPLED bulk formulations 
     148 
     149            ! albedo depends on cloud fraction because of non-linear spectral effects 
     150            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     151            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo  
     152            ! (zalb_ice) is computed within the bulk routine 
     153             
     154         END SELECT 
     155          
     156         !                                           ! Mask sea ice surface temperature 
    193157         DO jl = 1, jpl 
    194158            t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) ) 
    195159         END DO 
    196  
    197          IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) ) 
    198           
    199          IF( lk_cpl ) THEN 
    200             IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    201                ! 
    202                ! Compute mean albedo and temperature 
    203                zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) )  
    204                ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) )  
    205                ! 
    206             ENDIF 
    207          ENDIF 
    208                                                ! Bulk formulea - provides the following fields: 
     160      
     161         ! Bulk formulae  - provides the following fields: 
    209162         ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2] 
    210163         ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2] 
     
    215168         ! 
    216169         SELECT CASE( kblk ) 
    217          CASE( 3 )                                       ! CLIO bulk formulation 
    218             CALL blk_ice_clio( t_su , zalb_ice_cs, zalb_ice_os,                           & 
     170         CASE( jp_clio )                                       ! CLIO bulk formulation 
     171            CALL blk_ice_clio( t_su , zalb_cs    , zalb_os    , zalb_ice  ,               & 
    219172               &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   & 
    220173               &                      qla_ice    , dqns_ice   , dqla_ice  ,               & 
     
    222175               &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  ) 
    223176            !          
    224          CASE( 4 )                                       ! CORE bulk formulation 
    225             ! MV 2014 
    226             ! We must account for cloud fraction in the computation of the albedo 
    227             ! The present ref just uses the clear sky value 
    228             ! The overcast sky value is 0.06 higher, and polar skies are mostly overcast 
    229             ! CORE has no cloud fraction, hence we must prescribe it 
    230             ! Mean summer cloud fraction computed from CLIO = 0.81 
    231             zalb_ice(:,:,:) = 0.19 * zalb_ice_cs(:,:,:) + 0.81 * zalb_ice_os(:,:,:) 
    232             ! Following line, we replace zalb_ice_cs by simply zalb_ice 
     177            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     178               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     179 
     180         CASE( jp_core )                                       ! CORE bulk formulation 
    233181            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               & 
    234182               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   & 
     
    236184               &                      tprecip   , sprecip   ,                            & 
    237185               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  ) 
     186               ! 
     187            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     188               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
    238189            ! 
    239          CASE ( 5 ) 
    240             zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) ) 
     190         CASE ( jp_cpl ) 
    241191             
    242192            CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) 
    243193 
    244             CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=tn_ice    ) 
    245  
    246             ! Latent heat flux is forced to 0 in coupled : 
    247             !  it is included in qns (non-solar heat flux) 
    248             qla_ice  (:,:,:) = 0.0e0_wp 
    249             dqla_ice (:,:,:) = 0.0e0_wp 
     194            ! MV -> seb  
     195!           CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=t_su    ) 
     196 
     197!           IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     198!              &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     199!           ! Latent heat flux is forced to 0 in coupled : 
     200!           !  it is included in qns (non-solar heat flux) 
     201!           qla_ice  (:,:,:) = 0._wp 
     202!           dqla_ice (:,:,:) = 0._wp 
     203            ! END MV -> seb 
    250204            ! 
    251205         END SELECT 
    252  
    253          ! Average over all categories 
    254          IF( lk_cpl ) THEN 
    255          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    256  
    257             z_qns_ice_all  (:,:) = fice_ice_ave ( qns_ice  (:,:,:) ) 
    258             z_qsr_ice_all  (:,:) = fice_ice_ave ( qsr_ice  (:,:,:) ) 
    259             z_dqns_ice_all (:,:) = fice_ice_ave ( dqns_ice (:,:,:) ) 
    260             z_qla_ice_all  (:,:) = fice_ice_ave ( qla_ice  (:,:,:) ) 
    261             z_dqla_ice_all (:,:) = fice_ice_ave ( dqla_ice (:,:,:) ) 
    262  
    263             DO jl = 1, jpl 
    264                dqns_ice (:,:,jl) = z_dqns_ice_all (:,:) 
    265                dqla_ice (:,:,jl) = z_dqla_ice_all (:,:) 
    266             END DO 
    267             ! 
    268             IF ( ln_iceflx_ave ) THEN 
    269                DO jl = 1, jpl 
    270                   qns_ice  (:,:,jl) = z_qns_ice_all  (:,:) 
    271                   qsr_ice  (:,:,jl) = z_qsr_ice_all  (:,:) 
    272                   qla_ice  (:,:,jl) = z_qla_ice_all  (:,:) 
    273                END DO 
    274             END IF 
    275             ! 
    276             IF ( ln_iceflx_linear ) THEN 
    277                DO jl = 1, jpl 
    278                   qns_ice  (:,:,jl) = z_qns_ice_all(:,:) + z_dqns_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:)) 
    279                   qla_ice  (:,:,jl) = z_qla_ice_all(:,:) + z_dqla_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:)) 
    280                   qsr_ice  (:,:,jl) = (1.0e0_wp-zalb_ice(:,:,jl)) / (1.0e0_wp-zalb_ice_all(:,:)) * z_qsr_ice_all(:,:) 
    281                END DO 
    282             END IF 
    283          END IF 
    284          ENDIF 
     206          
    285207         !                                           !----------------------! 
    286208         !                                           ! LIM-3  time-stepping ! 
     
    290212         ! 
    291213         !                                           ! Store previous ice values 
    292 !!gm : remark   old_...   should becomes ...b  as tn versus tb   
    293          old_a_i  (:,:,:)   = a_i  (:,:,:)     ! ice area 
    294          old_e_i  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy 
    295          old_v_i  (:,:,:)   = v_i  (:,:,:)     ! ice volume 
    296          old_v_s  (:,:,:)   = v_s  (:,:,:)     ! snow volume  
    297          old_e_s  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy 
    298          old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content 
    299          old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content 
    300          old_u_ice(:,:)     = u_ice(:,:) 
    301          old_v_ice(:,:)     = v_ice(:,:) 
    302  
    303          ! trends    !!gm is it truly necessary ??? 
    304          d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp 
    305          d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp 
    306          d_e_i_thd  (:,:,:,:) = 0._wp   ;   d_e_i_trp  (:,:,:,:) = 0._wp 
    307          d_v_s_thd  (:,:,:)   = 0._wp   ;   d_v_s_trp  (:,:,:)   = 0._wp 
    308          d_e_s_thd  (:,:,:,:) = 0._wp   ;   d_e_s_trp  (:,:,:,:) = 0._wp 
    309          d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp 
    310          d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp 
    311          d_u_ice_dyn(:,:)     = 0._wp   ;   d_v_ice_dyn(:,:)     = 0._wp 
     214         a_i_b  (:,:,:)   = a_i  (:,:,:)     ! ice area 
     215         e_i_b  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy 
     216         v_i_b  (:,:,:)   = v_i  (:,:,:)     ! ice volume 
     217         v_s_b  (:,:,:)   = v_s  (:,:,:)     ! snow volume  
     218         e_s_b  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy 
     219         smv_i_b(:,:,:)   = smv_i(:,:,:)     ! salt content 
     220         oa_i_b (:,:,:)   = oa_i (:,:,:)     ! areal age content 
     221         u_ice_b(:,:)     = u_ice(:,:) 
     222         v_ice_b(:,:)     = v_ice(:,:) 
    312223 
    313224         ! salt, heat and mass fluxes 
    314225         sfx    (:,:) = 0._wp   ; 
    315          sfx_bri(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp  
     226         sfx_bri(:,:) = 0._wp   ;  
    316227         sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
    317228         sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
     
    334245         hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
    335246         hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
    336  
    337          ! 
    338          fhld  (:,:)    = 0._wp  
    339          fmmflx(:,:)    = 0._wp      
    340          ! part of solar radiation transmitted through the ice 
    341          ftr_ice(:,:,:) = 0._wp 
    342  
    343          ! diags 
    344          diag_trp_vi  (:,:) = 0._wp  ; diag_trp_vs(:,:) = 0._wp  ;  diag_trp_ei(:,:) = 0._wp  ;  diag_trp_es(:,:) = 0._wp 
    345          diag_heat_dhc(:,:) = 0._wp   
    346  
    347          ! dynamical invariants 
    348          delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp 
    349247 
    350248                          CALL lim_rst_opn( kt )     ! Open Ice restart file 
     
    372270         ENDIF 
    373271!                         !- Change old values for new values 
    374                           old_u_ice(:,:)   = u_ice (:,:) 
    375                           old_v_ice(:,:)   = v_ice (:,:) 
    376                           old_a_i(:,:,:)   = a_i (:,:,:) 
    377                           old_v_s(:,:,:)   = v_s (:,:,:) 
    378                           old_v_i(:,:,:)   = v_i (:,:,:) 
    379                           old_e_s(:,:,:,:) = e_s (:,:,:,:) 
    380                           old_e_i(:,:,:,:) = e_i (:,:,:,:) 
    381                           old_oa_i(:,:,:)  = oa_i(:,:,:) 
    382                           old_smv_i(:,:,:) = smv_i (:,:,:) 
     272                          u_ice_b(:,:)     = u_ice(:,:) 
     273                          v_ice_b(:,:)     = v_ice(:,:) 
     274                          a_i_b  (:,:,:)   = a_i (:,:,:) 
     275                          v_s_b  (:,:,:)   = v_s (:,:,:) 
     276                          v_i_b  (:,:,:)   = v_i (:,:,:) 
     277                          e_s_b  (:,:,:,:) = e_s (:,:,:,:) 
     278                          e_i_b  (:,:,:,:) = e_i (:,:,:,:) 
     279                          oa_i_b (:,:,:)   = oa_i (:,:,:) 
     280                          smv_i_b(:,:,:)   = smv_i(:,:,:) 
    383281  
    384282         ! ---------------------------------------------- 
     
    390288                          pfrld(:,:)   = 1._wp - at_i(:,:) 
    391289                          phicif(:,:)  = vt_i(:,:) 
     290 
     291                          ! MV -> seb 
     292                          SELECT CASE( kblk ) 
     293                             CASE ( jp_cpl ) 
     294                             CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su    ) 
     295                             IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
     296                          &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
     297                           ! Latent heat flux is forced to 0 in coupled : 
     298                           !  it is included in qns (non-solar heat flux) 
     299                             qla_ice  (:,:,:) = 0._wp 
     300                             dqla_ice (:,:,:) = 0._wp 
     301                          END SELECT 
     302                          ! END MV -> seb 
    392303                          ! 
    393304                          CALL lim_var_bv                 ! bulk brine volume (diag) 
     
    421332         IF( ln_nicep )   CALL lim_ctl( kt )              ! alerts in case of model crash 
    422333         ! 
     334         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice ) 
     335         ! 
    423336      ENDIF                                    ! End sea-ice time step only 
    424337 
     
    430343      !                                                   ! otherwise the atm.-ocean stresses are used everywhere 
    431344      IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents 
    432        
    433345!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    434       CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice ) 
    435  
    436       IF( lk_cpl ) THEN 
    437          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
    438             &    CALL wrk_dealloc( jpi, jpj, ztem_ice_all , zalb_ice_all , z_qsr_ice_all, z_qns_ice_all,   & 
    439             &                                z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
    440       ENDIF 
     346 
    441347      ! 
    442348      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim') 
    443349      ! 
    444350   END SUBROUTINE sbc_ice_lim 
    445  
    446  
     351    
     352    
     353      SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   & 
     354         &                          pdqn_ice, pqla_ice, pdql_ice, k_limflx ) 
     355      !!--------------------------------------------------------------------- 
     356      !!                  ***  ROUTINE sbc_ice_lim  *** 
     357      !!                    
     358      !! ** Purpose :   update the ice surface boundary condition by averaging and / or 
     359      !!                redistributing fluxes on ice categories                    
     360      !! 
     361      !! ** Method  :   average then redistribute  
     362      !! 
     363      !! ** Action  :    
     364      !!--------------------------------------------------------------------- 
     365      INTEGER                   , INTENT(in   ) ::   k_limflx   ! =-1 do nothing; =0 average ;  
     366                                                                ! =1 average and redistribute ; =2 redistribute 
     367      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature  
     368      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo 
     369      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqns_ice   ! non solar flux 
     370      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux 
     371      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity 
     372      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqla_ice   ! latent heat flux 
     373      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdql_ice   ! latent heat flux sensitivity 
     374      ! 
     375      INTEGER  ::   jl      ! dummy loop index 
     376      ! 
     377      REAL(wp), POINTER, DIMENSION(:,:) :: zalb_m    ! Mean albedo over all categories 
     378      REAL(wp), POINTER, DIMENSION(:,:) :: ztem_m    ! Mean temperature over all categories 
     379      ! 
     380      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m   ! Mean solar heat flux over all categories 
     381      REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m   ! Mean non solar heat flux over all categories 
     382      REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_m   ! Mean latent heat flux over all categories 
     383      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m   ! Mean d(qns)/dT over all categories 
     384      REAL(wp), POINTER, DIMENSION(:,:) :: z_dql_m   ! Mean d(qla)/dT over all categories 
     385      !!---------------------------------------------------------------------- 
     386 
     387      IF( nn_timing == 1 )  CALL timing_start('ice_lim_flx') 
     388      ! 
     389      ! 
     390      SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==! 
     391      CASE( 0 , 1 ) 
     392         CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m) 
     393         ! 
     394         z_qns_m(:,:) = fice_ice_ave ( pqns_ice (:,:,:) ) 
     395         z_qsr_m(:,:) = fice_ice_ave ( pqsr_ice (:,:,:) ) 
     396         z_dqn_m(:,:) = fice_ice_ave ( pdqn_ice (:,:,:) ) 
     397         z_qla_m(:,:) = fice_ice_ave ( pqla_ice (:,:,:) ) 
     398         z_dql_m(:,:) = fice_ice_ave ( pdql_ice (:,:,:) ) 
     399         DO jl = 1, jpl 
     400            pdqn_ice(:,:,jl) = z_dqn_m(:,:) 
     401            pdql_ice(:,:,jl) = z_dql_m(:,:) 
     402         END DO 
     403         ! 
     404         DO jl = 1, jpl 
     405            pqns_ice(:,:,jl) = z_qns_m(:,:) 
     406            pqsr_ice(:,:,jl) = z_qsr_m(:,:) 
     407            pqla_ice(:,:,jl) = z_qla_m(:,:) 
     408         END DO 
     409         ! 
     410         CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m) 
     411      END SELECT 
     412 
     413      SELECT CASE( k_limflx )                              !==  redistribution on all ice categories  ==! 
     414      CASE( 1 , 2 ) 
     415         CALL wrk_alloc( jpi,jpj, zalb_m, ztem_m ) 
     416         ! 
     417         zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) )  
     418         ztem_m(:,:) = fice_ice_ave ( ptn_ice  (:,:,:) )  
     419         DO jl = 1, jpl 
     420            pqns_ice(:,:,jl) = pqns_ice(:,:,jl) + pdqn_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:)) 
     421            pqla_ice(:,:,jl) = pqla_ice(:,:,jl) + pdql_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:)) 
     422            pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )  
     423         END DO 
     424         ! 
     425         CALL wrk_dealloc( jpi,jpj, zalb_m, ztem_m ) 
     426      END SELECT 
     427      ! 
     428      IF( nn_timing == 1 )  CALL timing_stop('ice_lim_flx') 
     429      ! 
     430   END SUBROUTINE ice_lim_flx 
     431    
     432    
    447433   SUBROUTINE lim_ctl( kt ) 
    448434      !!----------------------------------------------------------------------- 
     
    474460                  !WRITE(numout,*) ' at_i     ', at_i(ji,jj) 
    475461                  !WRITE(numout,*) ' Point - category', ji, jj, jl 
    476                   !WRITE(numout,*) ' a_i *** a_i_old ', a_i      (ji,jj,jl), old_a_i  (ji,jj,jl) 
    477                   !WRITE(numout,*) ' v_i *** v_i_old ', v_i      (ji,jj,jl), old_v_i  (ji,jj,jl) 
     462                  !WRITE(numout,*) ' a_i *** a_i_b   ', a_i      (ji,jj,jl), a_i_b  (ji,jj,jl) 
     463                  !WRITE(numout,*) ' v_i *** v_i_b   ', v_i      (ji,jj,jl), v_i_b  (ji,jj,jl) 
    478464                  !WRITE(numout,*) ' d_a_i_thd/trp   ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl) 
    479465                  !WRITE(numout,*) ' d_v_i_thd/trp   ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl) 
     
    585571               !DO jl = 1, jpl 
    586572                  !WRITE(numout,*) ' Category no: ', jl 
    587                   !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' old_a_i    : ', old_a_i  (ji,jj,jl)    
     573                  !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' a_i_b      : ', a_i_b  (ji,jj,jl)    
    588574                  !WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl) , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    589                   !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' old_v_i    : ', old_v_i  (ji,jj,jl)    
     575                  !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' v_i_b      : ', v_i_b  (ji,jj,jl)    
    590576                  !WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl) , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    591577                  !WRITE(numout,*) ' ' 
     
    676662      !!                n : number of the option 
    677663      !!------------------------------------------------------------------- 
    678       INTEGER         , INTENT(in) ::   kt      ! ocean time step 
     664      INTEGER         , INTENT(in) ::   kt            ! ocean time step 
    679665      INTEGER         , INTENT(in) ::   ki, kj, kn    ! ocean gridpoint indices 
    680666      CHARACTER(len=*), INTENT(in) ::   cd1           ! 
     
    763749               WRITE(numout,*) ' strength      : ', strength(ji,jj) 
    764750               WRITE(numout,*) ' d_u_ice_dyn   : ', d_u_ice_dyn(ji,jj), ' d_v_ice_dyn   : ', d_v_ice_dyn(ji,jj) 
    765                WRITE(numout,*) ' old_u_ice     : ', old_u_ice(ji,jj)  , ' old_v_ice     : ', old_v_ice(ji,jj)   
     751               WRITE(numout,*) ' u_ice_b       : ', u_ice_b(ji,jj)    , ' v_ice_b       : ', v_ice_b(ji,jj)   
    766752               WRITE(numout,*) 
    767753                
     
    773759                  WRITE(numout,*) ' t_su       : ', t_su(ji,jj,jl)             , ' t_s        : ', t_s(ji,jj,1,jl) 
    774760                  WRITE(numout,*) ' sm_i       : ', sm_i(ji,jj,jl)             , ' o_i        : ', o_i(ji,jj,jl) 
    775                   WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' old_a_i    : ', old_a_i(ji,jj,jl)    
     761                  WRITE(numout,*) ' a_i        : ', a_i(ji,jj,jl)              , ' a_i_b      : ', a_i_b(ji,jj,jl)    
    776762                  WRITE(numout,*) ' d_a_i_trp  : ', d_a_i_trp(ji,jj,jl)        , ' d_a_i_thd  : ', d_a_i_thd(ji,jj,jl)  
    777                   WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' old_v_i    : ', old_v_i(ji,jj,jl)    
     763                  WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' v_i_b      : ', v_i_b(ji,jj,jl)    
    778764                  WRITE(numout,*) ' d_v_i_trp  : ', d_v_i_trp(ji,jj,jl)        , ' d_v_i_thd  : ', d_v_i_thd(ji,jj,jl)  
    779                   WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' old_v_s    : ', old_v_s(ji,jj,jl)   
     765                  WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' v_s_b      : ', v_s_b(ji,jj,jl)   
    780766                  WRITE(numout,*) ' d_v_s_trp  : ', d_v_s_trp(ji,jj,jl)        , ' d_v_s_thd  : ', d_v_s_thd(ji,jj,jl) 
    781                   WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' old_ei1    : ', old_e_i(ji,jj,1,jl)/1.0e9  
     767                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' ei1        : ', e_i_b(ji,jj,1,jl)/1.0e9  
    782768                  WRITE(numout,*) ' de_i1_trp  : ', d_e_i_trp(ji,jj,1,jl)/1.0e9, ' de_i1_thd  : ', d_e_i_thd(ji,jj,1,jl)/1.0e9 
    783                   WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' old_ei2    : ', old_e_i(ji,jj,2,jl)/1.0e9   
     769                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' ei2_b      : ', e_i_b(ji,jj,2,jl)/1.0e9   
    784770                  WRITE(numout,*) ' de_i2_trp  : ', d_e_i_trp(ji,jj,2,jl)/1.0e9, ' de_i2_thd  : ', d_e_i_thd(ji,jj,2,jl)/1.0e9 
    785                   WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' old_e_snow : ', old_e_s(ji,jj,1,jl)  
     771                  WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl)  
    786772                  WRITE(numout,*) ' d_e_s_trp  : ', d_e_s_trp(ji,jj,1,jl)      , ' d_e_s_thd  : ', d_e_s_thd(ji,jj,1,jl) 
    787                   WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' old_smv_i  : ', old_smv_i(ji,jj,jl)    
     773                  WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' smv_i_b    : ', smv_i_b(ji,jj,jl)    
    788774                  WRITE(numout,*) ' d_smv_i_trp: ', d_smv_i_trp(ji,jj,jl)      , ' d_smv_i_thd: ', d_smv_i_thd(ji,jj,jl)  
    789                   WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' old_oa_i   : ', old_oa_i(ji,jj,jl) 
     775                  WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' oa_i_b     : ', oa_i_b(ji,jj,jl) 
    790776                  WRITE(numout,*) ' d_oa_i_trp : ', d_oa_i_trp(ji,jj,jl)       , ' d_oa_i_thd : ', d_oa_i_thd(ji,jj,jl) 
    791777               END DO !jl 
     
    795781               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    796782               WRITE(numout,*) ' - Heat fluxes in and out the ice ***' 
    797                WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( old_a_i(ji,jj,:) * qsr_ice(ji,jj,:) ) 
    798                WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( old_a_i(ji,jj,:) * qns_ice(ji,jj,:) ) 
     783               WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) ) 
     784               WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) ) 
    799785               WRITE(numout,*) 
    800786               WRITE(numout,*)  
     
    854840         END DO 
    855841      END DO 
    856  
     842      ! 
    857843   END SUBROUTINE lim_prt_state 
     844    
     845      
     846   FUNCTION fice_cell_ave ( ptab ) 
     847      !!-------------------------------------------------------------------------- 
     848      !! * Compute average over categories, for grid cell (ice covered and free ocean) 
     849      !!-------------------------------------------------------------------------- 
     850      REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave 
     851      REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab 
     852      INTEGER :: jl ! Dummy loop index 
     853       
     854      fice_cell_ave (:,:) = 0.0_wp 
     855       
     856      DO jl = 1, jpl 
     857         fice_cell_ave (:,:) = fice_cell_ave (:,:) & 
     858            &                  + a_i (:,:,jl) * ptab (:,:,jl) 
     859      END DO 
     860       
     861   END FUNCTION fice_cell_ave 
     862    
     863    
     864   FUNCTION fice_ice_ave ( ptab ) 
     865      !!-------------------------------------------------------------------------- 
     866      !! * Compute average over categories, for ice covered part of grid cell 
     867      !!-------------------------------------------------------------------------- 
     868      REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave 
     869      REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab 
     870 
     871      fice_ice_ave (:,:) = 0.0_wp 
     872      WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:) 
     873 
     874   END FUNCTION fice_ice_ave 
     875 
    858876 
    859877#else 
Note: See TracChangeset for help on using the changeset viewer.