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

Ignore:
Timestamp:
2014-11-27T16:28:53+01:00 (9 years ago)
Author:
cetlod
Message:

2014/dev_CNRS_2014 : Merge in the trunk changes between 4674 and 4728, see ticket #1415

File:
1 edited

Legend:

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

    r4896 r4900  
    5959   USE prtctl          ! Print control 
    6060   USE lib_fortran     !  
     61   USE cpl_oasis3, ONLY : lk_cpl 
    6162 
    6263#if defined key_bdy  
     
    6869 
    6970   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90 
     71   PUBLIC lim_prt_state 
    7072    
    7173   !! * Substitutions 
     
    133135      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE) 
    134136      !! 
    135       INTEGER  ::   jl      ! dummy loop index 
     137      INTEGER  ::   ji, jj, jl, jk      ! dummy loop index 
    136138      REAL(wp) ::   zcoef   ! local scalar 
    137139      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky 
     
    146148      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories 
    147149      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   ! 
    148152      !!---------------------------------------------------------------------- 
    149153 
     
    152156      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
    153157 
    154       CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
    155  
    156 #if defined key_coupled 
    157       IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_alloc( jpi,jpj,jpl, zalb_ice) 
    158       IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
    159          &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
    160 #endif 
     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, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
     163      ENDIF 
    161164 
    162165      IF( kt == nit000 ) THEN 
     
    168171         ! 
    169172         IF( ln_nicep ) THEN      ! control print at a given point 
    170             jiindx = 177   ;   jjindx = 112 
     173            jiindx = 15    ;   jjindx =  44 
    171174            IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    172175         ENDIF 
     
    176179      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  ! 
    177180         !                                     !----------------------! 
    178          !                                           !  Bulk Formulea ! 
     181         !                                           !  Bulk Formulae ! 
    179182         !                                           !----------------! 
    180183         ! 
     
    182185         v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean) 
    183186         ! 
    184          t_bo(:,:) = eos_fzp( sss_m ) +  rt0         ! masked sea surface freezing temperature [Kelvin] 
    185          !                                           ! (set to rt0 over land) 
     187         t_bo(:,:) = ( eos_fzp( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) )   ! masked sea surface freezing temperature [Kelvin] 
     188         !                                                                                      ! (set to rt0 over land) 
    186189         CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo 
    187190 
     
    192195         IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) ) 
    193196          
    194 #if defined key_coupled 
    195          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    196             ! 
    197             ! Compute mean albedo and temperature 
    198             zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) )  
    199             ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) )  
    200             ! 
     197         IF( lk_cpl ) THEN 
     198            IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
     199               ! 
     200               ! Compute mean albedo and temperature 
     201               zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) )  
     202               ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) )  
     203               ! 
     204            ENDIF 
    201205         ENDIF 
    202 #endif 
    203206                                               ! Bulk formulea - provides the following fields: 
    204207         ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2] 
     
    218221            !          
    219222         CASE( 4 )                                       ! CORE bulk formulation 
    220             CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice_cs,               & 
     223            ! MV 2014 
     224            ! We must account for cloud fraction in the computation of the albedo 
     225            ! The present ref just uses the clear sky value 
     226            ! The overcast sky value is 0.06 higher, and polar skies are mostly overcast 
     227            ! CORE has no cloud fraction, hence we must prescribe it 
     228            ! Mean summer cloud fraction computed from CLIO = 0.81 
     229            zalb_ice(:,:,:) = 0.19 * zalb_ice_cs(:,:,:) + 0.81 * zalb_ice_os(:,:,:) 
     230            ! Following line, we replace zalb_ice_cs by simply zalb_ice 
     231            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               & 
    221232               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   & 
    222233               &                      qla_ice   , dqns_ice  , dqla_ice   ,               & 
     
    239250 
    240251         ! Average over all categories 
    241 #if defined key_coupled 
     252         IF( lk_cpl ) THEN 
    242253         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    243254 
     
    269280            END IF 
    270281         END IF 
    271 #endif 
     282         ENDIF 
    272283         !                                           !----------------------! 
    273284         !                                           ! LIM-3  time-stepping ! 
     
    285296         old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content 
    286297         old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content 
    287          ! 
    288          old_u_ice(:,:) = u_ice(:,:) 
    289          old_v_ice(:,:) = v_ice(:,:) 
    290          !                                           ! intialisation to zero    !!gm is it truly necessary ??? 
     298         old_u_ice(:,:)     = u_ice(:,:) 
     299         old_v_ice(:,:)     = v_ice(:,:) 
     300 
     301         ! trends    !!gm is it truly necessary ??? 
    291302         d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp 
    292303         d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp 
     
    296307         d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp 
    297308         d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp 
    298          ! 
    299          d_u_ice_dyn(:,:) = 0._wp 
    300          d_v_ice_dyn(:,:) = 0._wp 
    301          ! 
    302          sfx    (:,:) = 0._wp   ;   sfx_thd  (:,:) = 0._wp 
    303          sfx_bri(:,:) = 0._wp   ;   sfx_mec  (:,:) = 0._wp   ;   sfx_res  (:,:) = 0._wp 
    304          fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp 
    305          fhmec  (:,:) = 0._wp   ;    
    306          fmmec  (:,:) = 0._wp 
    307          fmmflx (:,:) = 0._wp      
    308          focea2D(:,:) = 0._wp 
    309          fsup2D (:,:) = 0._wp 
    310  
    311          ! used in limthd.F90 
    312          rdvosif(:,:) = 0._wp   ! variation of ice volume at surface 
    313          rdvobif(:,:) = 0._wp   ! variation of ice volume at bottom 
    314          fdvolif(:,:) = 0._wp   ! total variation of ice volume 
    315          rdvonif(:,:) = 0._wp   ! lateral variation of ice volume 
    316          fstric (:,:) = 0._wp   ! part of solar radiation transmitted through the ice 
    317          ffltbif(:,:) = 0._wp   ! linked with fstric 
    318          qfvbq  (:,:) = 0._wp   ! linked with fstric 
    319          rdm_snw(:,:) = 0._wp   ! variation of snow mass per unit area 
    320          rdm_ice(:,:) = 0._wp   ! variation of ice mass per unit area 
    321          hicifp (:,:) = 0._wp   ! daily thermodynamic ice production.  
    322          ! 
    323          diag_sni_gr(:,:) = 0._wp   ;   diag_lat_gr(:,:) = 0._wp 
    324          diag_bot_gr(:,:) = 0._wp   ;   diag_dyn_gr(:,:) = 0._wp 
    325          diag_bot_me(:,:) = 0._wp   ;   diag_sur_me(:,:) = 0._wp 
    326          diag_res_pr(:,:) = 0._wp   ;   diag_trp_vi(:,:) = 0._wp 
     309         d_u_ice_dyn(:,:)     = 0._wp   ;   d_v_ice_dyn(:,:)     = 0._wp 
     310 
     311         ! salt, heat and mass fluxes 
     312         sfx    (:,:) = 0._wp   ; 
     313         sfx_bri(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp  
     314         sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
     315         sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
     316         sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
     317         sfx_res(:,:) = 0._wp 
     318 
     319         wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
     320         wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
     321         wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
     322         wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
     323         wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
     324         wfx_spr(:,:) = 0._wp   ;    
     325 
     326         hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
     327         hfx_thd(:,:) = 0._wp   ;    
     328         hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
     329         hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
     330         hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
     331         hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
     332         hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
     333         hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
     334 
     335         ! 
     336         fhld  (:,:)    = 0._wp  
     337         fmmflx(:,:)    = 0._wp      
     338         ! part of solar radiation transmitted through the ice 
     339         ftr_ice(:,:,:) = 0._wp 
     340 
     341         ! diags 
     342         diag_trp_vi  (:,:) = 0._wp  ; diag_trp_vs(:,:) = 0._wp  ;  diag_trp_ei(:,:) = 0._wp  ;  diag_trp_es(:,:) = 0._wp 
     343         diag_heat_dhc(:,:) = 0._wp   
     344 
    327345         ! dynamical invariants 
    328346         delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp 
     
    375393                          zcoef = rdt_ice /rday           !  Ice natural aging 
    376394                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 
    377                           CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin) 
    378395         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
    379396                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  ! 
     
    391408         !                                           ! Diagnostics and outputs  
    392409         IF (ln_limdiaout) CALL lim_diahsb 
    393 !clem # if ! defined key_iomput 
     410 
    394411                          CALL lim_wri( 1  )              ! Ice outputs  
    395 !clem # endif 
     412 
    396413         IF( kt == nit000 .AND. ln_rstart )   & 
    397414            &             CALL iom_close( numrir )        ! clem: close input ice restart file 
     
    413430       
    414431!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    415       ! 
    416       CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
    417  
    418 #if defined key_coupled 
    419       IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice) 
    420       IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
    421          &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
    422 #endif 
     432      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice ) 
     433 
     434      IF( lk_cpl ) THEN 
     435         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
     436            &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
     437      ENDIF 
    423438      ! 
    424439      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim') 
     
    534549!                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    535550!                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    536 !                 WRITE(numout,*) ' s_i_newice           : ', s_i_newice(ji,jj,1:jpl) 
    537551!                 WRITE(numout,*)  
    538552                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     
    591605               !WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
    592606               !WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
    593                !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) 
    594                !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) 
    595                !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice 
    596                !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice 
    597                !WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj) 
    598                !WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj) 
    599                !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice 
    600                !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice 
    601                !WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)  
    602                !WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj)  
    603                !WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj)  
    604                !WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)  
    605                !WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj)  
    606607               ! 
    607608               !CALL lim_prt_state( kt, ji, jj, 2, '   ') 
     
    790791               WRITE(numout,*) ' - Heat / FW fluxes ' 
    791792               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    792                WRITE(numout,*) ' emp        : ', emp      (ji,jj) 
    793                WRITE(numout,*) ' sfx        : ', sfx      (ji,jj) 
    794                WRITE(numout,*) ' sfx_thd    : ', sfx_thd(ji,jj) 
    795                WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ji,jj) 
    796                WRITE(numout,*) ' sfx_mec    : ', sfx_mec  (ji,jj) 
    797                WRITE(numout,*) ' sfx_res    : ', sfx_res(ji,jj) 
    798                WRITE(numout,*) ' fmmec      : ', fmmec    (ji,jj) 
    799                WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj) 
    800                WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj) 
    801                WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj) 
     793               WRITE(numout,*) ' - Heat fluxes in and out the ice ***' 
     794               WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( old_a_i(ji,jj,:) * qsr_ice(ji,jj,:) ) 
     795               WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( old_a_i(ji,jj,:) * qns_ice(ji,jj,:) ) 
     796               WRITE(numout,*) 
    802797               WRITE(numout,*)  
    803798               WRITE(numout,*) ' sst        : ', sst_m(ji,jj)   
     
    829824               WRITE(numout,*) ' qsr       : ', qsr(ji,jj) 
    830825               WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    831                WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
    832                WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) * r1_rdtice 
    833                WRITE(numout,*) ' qldif     : ', qldif(ji,jj) * r1_rdtice 
     826               WRITE(numout,*) 
     827               WRITE(numout,*) ' hfx_mass     : ', hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_snw(ji,jj) + hfx_res(ji,jj) 
     828               WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj) 
     829               WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj) 
     830               WRITE(numout,*) ' dhc          : ', diag_heat_dhc(ji,jj)               
     831               WRITE(numout,*) 
     832               WRITE(numout,*) ' hfx_dyn      : ', hfx_dyn(ji,jj) 
     833               WRITE(numout,*) ' hfx_thd      : ', hfx_thd(ji,jj) 
     834               WRITE(numout,*) ' hfx_res      : ', hfx_res(ji,jj) 
     835               WRITE(numout,*) ' fhtur        : ', fhtur(ji,jj)  
     836               WRITE(numout,*) ' qlead        : ', qlead(ji,jj) * r1_rdtice 
    834837               WRITE(numout,*) 
    835838               WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
    836839               WRITE(numout,*) ' emp       : ', emp    (ji,jj) 
    837                WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
    838840               WRITE(numout,*) ' sfx       : ', sfx    (ji,jj) 
    839841               WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj) 
    840                WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ji,jj) 
    841                WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
    842                WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
     842               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
     843               WRITE(numout,*) ' sfx_dyn   : ', sfx_dyn(ji,jj) 
    843844               WRITE(numout,*) 
    844845               WRITE(numout,*) ' - Momentum fluxes ' 
    845846               WRITE(numout,*) ' utau      : ', utau(ji,jj)  
    846847               WRITE(numout,*) ' vtau      : ', vtau(ji,jj) 
    847             ENDIF 
     848            ENDIF  
    848849            WRITE(numout,*) ' ' 
    849850            ! 
Note: See TracChangeset for help on using the changeset viewer.