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 13886 for NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icethd.F90 – NEMO

Ignore:
Timestamp:
2020-11-26T15:24:38+01:00 (3 years ago)
Author:
rlod
Message:

phasing with trunk at revision r13787

Location:
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source

    • Property svn:externals
      •  

        old new  
        88 
        99# SETTE 
        10 ^/utils/CI/sette@13507        sette 
         10^/utils/CI/sette@13559        sette 
  • NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icethd.F90

    r13472 r13886  
    1818   USE ice            ! sea-ice: variables 
    1919!!gm list trop longue ==>>> why not passage en argument d'appel ? 
    20    USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl 
     20   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 
    2121   USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 
    2222      &                 qml_ice, qcn_ice, qtr_ice_top 
     
    3030   USE icethd_pnd     ! sea-ice: melt ponds 
    3131   USE iceitd         ! sea-ice: remapping thickness distribution 
     32   USE icecor         ! sea-ice: corrections 
    3233   USE icetab         ! sea-ice: 1D <==> 2D transformation 
    3334   USE icevar         ! sea-ice: operations 
     
    5253   LOGICAL ::   ln_icedO         ! activate ice growth in open-water (T) or not (F) 
    5354   LOGICAL ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F) 
    54    LOGICAL ::   ln_leadhfx       !  heat in the leads is used to melt sea-ice before warming the ocean 
     55   LOGICAL ::   ln_leadhfx       ! heat in the leads is used to melt sea-ice before warming the ocean 
    5556 
    5657   !! for convergence tests 
     
    9192      ! 
    9293      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices 
    93       REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg 
    94       REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04) 
    95       REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient 
    96       REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 
     94      REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos 
     95      REAL(wp), PARAMETER :: zfric_umin = 0._wp       ! lower bound for the friction velocity (cice value=5.e-04) 
     96      REAL(wp), PARAMETER :: zch        = 0.0057_wp   ! heat transfer coefficient 
     97      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric, zvel   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 
    9798      ! 
    9899      !!------------------------------------------------------------------- 
     
    124125               &                    (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
    125126               &                     + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 
     127            zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + & 
     128               &                         ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 
    126129         END_2D 
    127130      ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean 
     
    130133               &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
    131134               &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 
     135            zvel(ji,jj) = 0._wp 
    132136         END_2D 
    133137      ENDIF 
    134       CALL lbc_lnk( 'icethd', zfric, 'T', 1.0_wp ) 
     138      CALL lbc_lnk_multi( 'icethd', zfric, 'T',  1.0_wp, zvel, 'T', 1.0_wp ) 
    135139      ! 
    136140      !--------------------------------------------------------------------! 
     
    140144         rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 
    141145         ! 
    142          !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
    143          !           !  practically no "direct lateral ablation" 
    144          !            
    145          !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean  
    146          !           !  temperature and turbulent mixing (McPhee, 1992) 
    147          ! 
    148146         ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 
    149147         zqld =  tmask(ji,jj,1) * rDt_ice *  & 
     
    151149            &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 
    152150 
    153          ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- ! 
     151         ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 
     152         !     (mostly<0 but >0 if supercooling) 
    154153         zqfr     = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1)  ! both < 0 (t_bo < sst) and > 0 (t_bo > sst) 
    155154         zqfr_neg = MIN( zqfr , 0._wp )                                                                    ! only < 0 
    156  
    157          ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2) 
     155         zqfr_pos = MAX( zqfr , 0._wp )                                                                    ! only > 0 
     156 
     157         ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 
     158         !     (mostly>0 but <0 if supercooling) 
    158159         zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin )  
    159          qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 
    160  
    161          qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 
     160         qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 
     161          
    162162         ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach  
    163163         !                              the freezing point, so that we do not have SST < T_freeze 
    164          !                              This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0 
    165  
    166          !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice 
    167          qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 
    168  
    169          ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting  
    170          ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget 
    171          IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN 
    172             IF( ln_leadhfx ) THEN   ;   fhld(ji,jj) = rswitch * zqld * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
    173             ELSE                    ;   fhld(ji,jj) = 0._wp 
    174             ENDIF 
     164         !                              This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 
     165         !                              The following formulation is ok for both normal conditions and supercooling 
     166         qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 
     167 
     168         ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 
     169         !     qlead is the energy received from the atm. in the leads. 
     170         !     If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld  (W/m2) 
     171         !     If cooling (zqld <  0), then the energy in the leads is used to grow ice in open water    => qlead (J.m-2) 
     172         IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
     173            ! upper bound for fhld: fhld should be equal to zqld 
     174            !                        but we have to make sure that this heat will not make the sst drop below the freezing point 
     175            !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 
     176            !                        The following formulation is ok for both normal conditions and supercooling 
     177            fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) &  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
     178               &                                 - qsb_ice_bot(ji,jj) ) 
    175179            qlead(ji,jj) = 0._wp 
    176180         ELSE 
    177181            fhld (ji,jj) = 0._wp 
     182            ! upper bound for qlead: qlead should be equal to zqld 
     183            !                        but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 
     184            !                        The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 
     185            !                        and freezing point is reached if zqfr = zqld - qsb*a/dt 
     186            !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 
     187            !                        The following formulation is ok for both normal conditions and supercooling 
     188            qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 
    178189         ENDIF 
    179190         ! 
    180          ! Net heat flux on top of the ice-ocean [W.m-2] 
    181          ! --------------------------------------------- 
    182          qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj)  
     191         ! If ice is landfast and ice concentration reaches its max 
     192         ! => stop ice formation in open water 
     193         IF(  zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 )   qlead(ji,jj) = 0._wp 
     194         ! 
     195         ! If the grid cell is almost fully covered by ice (no leads) 
     196         ! => stop ice formation in open water 
     197         IF( at_i(ji,jj) >= (1._wp - epsi10) )   qlead(ji,jj) = 0._wp 
     198         ! 
     199         ! If ln_leadhfx is false 
     200         ! => do not use energy of the leads to melt sea-ice 
     201         IF( .NOT.ln_leadhfx )   fhld(ji,jj) = 0._wp 
     202         ! 
    183203      END_2D 
    184204       
     
    191211      ENDIF 
    192212 
    193       ! --------------------------------------------------------------------- 
    194       ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2] 
    195       ! --------------------------------------------------------------------- 
    196       !     First  step here              :  non solar + precip - qlead - qsensible 
    197       !     Second step in icethd_dh      :  heat remaining if total melt (zq_rema)  
    198       !     Third  step in iceupdate.F90  :  heat from ice-ocean mass exchange (zf_mass) + solar 
    199       qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean                
    200          &             - qlead(:,:) * r1_Dt_ice                                &  ! heat flux taken from the ocean where there is open water ice formation 
    201          &             - at_i (:,:) * qsb_ice_bot(:,:)                         &  ! heat flux taken by sensible flux 
    202          &             - at_i (:,:) * fhld       (:,:)                            ! heat flux taken during bottom growth/melt  
    203       !                                                                           !    (fhld should be 0 while bott growth) 
    204213      !-------------------------------------------------------------------------------------------! 
    205214      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 
     
    231240                              CALL ice_thd_dh                           ! Ice-Snow thickness    
    232241                              CALL ice_thd_pnd                          ! Melt ponds formation 
    233                               CALL ice_thd_ent( e_i_1d(1:npti,:), .true. )      ! Ice enthalpy remapping 
     242                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping 
    234243            ENDIF 
    235244                              CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !     
     
    254263      ! 
    255264      IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- ! 
     265      ! 
     266                              CALL ice_cor( kt , 2 )                ! --- Corrections --- ! 
     267      ! 
     268      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice              ! ice natural aging incrementation      
    256269      ! 
    257270      ! convergence tests 
     
    418431         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res       ) 
    419432         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif   ) 
    420          CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai     ) 
    421433         ! 
    422434         ! ocean surface fields 
    423435         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m ) 
    424436         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m ) 
     437         CALL tab_2d_1d( npti, nptidx(1:npti), frq_m_1d(1:npti), frq_m ) 
    425438         ! 
    426439         ! to update ice age 
     
    510523         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res     ) 
    511524         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif ) 
    512          CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai   ) 
    513525         ! 
    514526         CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d    (1:npti), qns_ice    (:,:,kl) ) 
Note: See TracChangeset for help on using the changeset viewer.