Changeset 8623


Ignore:
Timestamp:
2017-10-12T23:24:57+02:00 (4 years ago)
Author:
clem
Message:

rewrite melt ponds

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90

    r8565 r8623  
    8989 
    9090   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip 
    91    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_1d        !: <==> the 2D  at_i 
     91   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_1d       !: <==> the 2D  at_i 
    9292   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ato_i_1d      !: <==> the 2D  ato_i 
    9393   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhtur_1d      !: <==> the 2D  fhtur 
     
    103103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_i_1d        !: <==> the 2D  a_i 
    104104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ib_1d       !: <==> the 2D  a_i_b 
    105    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_i_1d       !: 
    106    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ib_1d      !: 
    107    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_s_1d       !: 
     105   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_i_1d        !: 
     106   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ib_1d       !: 
     107   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_s_1d        !: 
    108108   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_su         !: Surface Conduction flux  
    109109   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_bo_i       !: Bottom  Conduction flux  
     
    112112   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_sub      !: Ice surface sublimation [m] 
    113113   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott     !: Ice bottom accretion/ablation  [m] 
     114   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_mlt      !: Snow melt [m] 
    114115   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice    !: Snow ice formation             [m of ice] 
    115    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_1d       !: Ice bulk salinity [ppt] 
     116   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_1d        !: Ice bulk salinity [ppt] 
    116117   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new       !: Salinity of new ice at the bottom 
     118   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_i_1d        !: 
     119   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_s_1d        !: 
     120   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sv_i_1d       !: 
     121   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   oa_i_1d       !: 
    117122   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d       !: 
    118123   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_ip_1d       !: 
    119    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_i_1d        !: 
    120    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_s_1d        !: 
    121    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sv_i_1d      !: 
    122    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   oa_i_1d       !: 
    123  
    124    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d   !: corresponding to the 2D var  t_s 
    125    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_i_1d   !: corresponding to the 2D var  t_i 
    126    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sz_i_1d  !: profiled ice salinity 
    127    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e_i_1d   !:    Ice  enthalpy per unit volume 
    128    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e_s_1d   !:    Snow enthalpy per unit volume 
    129  
    130    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   eh_i_old !: ice heat content (q*h, J.m-2) 
    131    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_old  !: ice thickness layer (m) 
     124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ip_1d       !: 
     125   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !: 
     126 
     127   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s 
     128   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_i_1d      !: corresponding to the 2D var  t_i 
     129   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sz_i_1d     !: profiled ice salinity 
     130   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e_i_1d      !:    Ice  enthalpy per unit volume 
     131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e_s_1d      !:    Snow enthalpy per unit volume 
     132 
     133   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   eh_i_old    !: ice heat content (q*h, J.m-2) 
     134   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_old     !: ice thickness layer (m) 
    132135 
    133136   ! Conduction flux diagnostics (SIMIP) 
     
    200203         &      h_i_1d  (jpij) , h_ib_1d  (jpij) , h_s_1d (jpij) , fc_su  (jpij) , fc_bo_i(jpij) ,  &     
    201204         &      dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) ,                                  &     
    202          &      dh_i_bott(jpij) , dh_snowice(jpij) , s_i_1d (jpij) , s_i_new(jpij) , & 
     205         &      dh_i_bott(jpij) , dh_s_mlt(jpij), dh_snowice(jpij) , s_i_1d (jpij) , s_i_new(jpij) , & 
    203206         &      a_ip_1d  (jpij) , v_ip_1d   (jpij) , v_i_1d  (jpij) , v_s_1d (jpij) , & 
     207         &      h_ip_1d  (jpij) , a_ip_frac_1d(jpij) , & 
    204208         &      sv_i_1d (jpij) , oa_i_1d   (jpij) , STAT=ierr(ii) ) 
    205209      ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90

    r8567 r8623  
    7676      IF( icount == 0 ) THEN 
    7777         !                          ! water flux 
    78          pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
    79             &                    wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
    80             &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    & 
     78         pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  & 
     79            &                    wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  & 
     80            &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  & 
     81            &                    wfx_ice_sub(:,:) + wfx_spr(:,:)  & 
    8182            &                  ) * e1e2t(:,:) ) * zconv 
    8283         ! 
     
    9697 
    9798         pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
    98             &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv 
     99            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv 
    99100 
    100101      ELSEIF( icount == 1 ) THEN 
    101102 
    102103         ! water flux 
    103          zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
    104             &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
    105             &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        & 
     104         zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  & 
     105            &                wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  & 
     106            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  & 
     107            &                wfx_ice_sub(:,:) + wfx_spr(:,:)  & 
    106108            &              ) * e1e2t(:,:) ) * zconv - pdiag_fv 
    107109 
     
    117119  
    118120         ! outputs 
    119          zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  & 
     121         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv  & 
    120122            &     - pdiag_v ) * r1_rdtice - zfv ) * rday 
    121123 
    122          zs = ( ( glob_sum( SUM( sv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  & 
     124         zs = ( ( glob_sum( SUM( sv_i * rhoic             , dim=3 ) * e1e2t ) * zconv  & 
    123125            &     - pdiag_s ) * r1_rdtice + zfs ) * rday 
    124126 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

    r8598 r8623  
    231231            s_i_new   (1:npti) = 0._wp ; dh_s_tot (1:npti) = 0._wp  ! --- some init --- !  (important to have them here)  
    232232            dh_i_surf (1:npti) = 0._wp ; dh_i_bott(1:npti) = 0._wp 
    233             dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp 
     233            dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 
    234234            ! 
    235235            IF( ln_icedH ) THEN                                     ! --- growing/melting --- ! 
    236236                              CALL ice_thd_zdf                             ! Ice/Snow Temperature profile 
    237237                              CALL ice_thd_dh                              ! Ice/Snow thickness    
     238                              CALL ice_thd_pnd                             ! Melt ponds formation 
    238239                              CALL ice_thd_ent( e_i_1d(1:npti,:) )         ! Ice enthalpy remapping 
    239240            ENDIF 
     
    273274      ! 
    274275      IF( ln_icedO )       CALL ice_thd_do                 ! --- frazil ice growing in leads --- ! 
    275       ! 
    276                            CALL ice_thd_pnd( kt )          ! --- Melt ponds --- ! 
    277276      ! 
    278277      ! controls 
     
    365364         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             ) 
    366365         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     ) 
    367          CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl)     ) 
    368          CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl)     ) 
     366         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     ) 
     367         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     ) 
    369368         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     ) 
    370          CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl)     ) 
     369         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     ) 
    371370         DO jk = 1, nlay_s 
    372             CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)   ) 
    373             CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)   ) 
     371            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    ) 
     372            CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    ) 
    374373         END DO 
    375374         DO jk = 1, nlay_i 
    376             CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl)   ) 
    377             CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl)   ) 
    378             CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)   ) 
    379          END DO 
     375            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  ) 
     376            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  ) 
     377            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  ) 
     378         END DO 
     379         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) ) 
     380         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
     381         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 
    380382         ! 
    381383         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice        ) 
     
    406408         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr          ) 
    407409         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam          ) 
     410         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          ) 
    408411         ! 
    409412         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog          ) 
     
    457460         ! 
    458461         ! Change thickness to volume (replaces routine ice_var_eqv2glo) 
    459          v_i_1d(1:npti)  = h_i_1d(1:npti) * a_i_1d(1:npti) 
    460          v_s_1d(1:npti)  = h_s_1d(1:npti) * a_i_1d(1:npti) 
    461          sv_i_1d(1:npti) = s_i_1d(1:npti) * v_i_1d(1:npti) 
     462         v_i_1d (1:npti) = h_i_1d (1:npti) * a_i_1d (1:npti) 
     463         v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti) 
     464         sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti) 
     465         v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 
    462466          
    463467         CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             ) 
    464468         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     ) 
    465          CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl)     ) 
    466          CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl)     ) 
     469         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     ) 
     470         CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     ) 
    467471         CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     ) 
    468          CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl)     ) 
     472         CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     ) 
    469473         DO jk = 1, nlay_s 
    470             CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) ) 
    471             CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) ) 
     474            CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    ) 
     475            CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    ) 
    472476         END DO 
    473477         DO jk = 1, nlay_i 
    474             CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl) ) 
    475             CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl) ) 
    476             CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 
    477          END DO 
     478            CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  ) 
     479            CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  ) 
     480            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  ) 
     481         END DO 
     482         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) ) 
     483         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
     484         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 
    478485         ! 
    479486         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) 
     
    491498         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        ) 
    492499         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        ) 
     500         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
    493501         ! 
    494502         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        ) 
     
    526534         CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) ) 
    527535         CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) ) 
     536         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) ) 
    528537         ! 
    529538      END SELECT 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

    r8565 r8623  
    148148            ! Contribution to mass flux 
    149149            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * h_s_1d(ji) * a_i_1d(ji) * r1_rdtice 
     150            dh_s_mlt(ji)       = dh_s_mlt(ji) - h_s_1d(ji) 
    150151            ! updates 
    151             h_s_1d(ji)   = 0._wp 
     152            h_s_1d(ji)    = 0._wp 
    152153            e_s_1d (ji,1) = 0._wp 
    153154            t_s_1d (ji,1) = rt0 
     
    211212         ! snow melting only = water into the ocean (then without snow precip), >0 
    212213         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
     214         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    213215         ! updates available heat + precipitations after melting 
    214216         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )       
     
    233235            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice  
    234236            ! snow melting only = water into the ocean (then without snow precip) 
    235             wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     237            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     238            dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,jk) 
    236239            ! updates available heat + thickness 
    237             zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
     240            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
    238241            h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 
    239242         END DO 
     
    571574         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 
    572575         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1) 
    573          h_s_1d   (ji)  = h_s_1d(ji)   + zdeltah(ji,1) 
     576         h_s_1d   (ji)   = h_s_1d(ji)   + zdeltah(ji,1) 
    574577         
    575578         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2) 
     
    577580         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 
    578581         ! Contribution to mass flux 
    579          wfx_snw_sum_1d(ji)  =  wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
     582         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
     583         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    580584         !     
    581585         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90

    r8598 r8623  
    66   !! history :       ! Original code by Daniela Flocco and Adrian Turner 
    77   !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1) 
    8    !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6 
     8   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO4 
    99   !!---------------------------------------------------------------------- 
    1010#if defined key_lim3 
    1111   !!---------------------------------------------------------------------- 
    12    !!   'key_lim3' :                                 LIM3 sea-ice model 
     12   !!   'key_lim3' :                                     ESIM sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    1414   !!   ice_thd_pnd_init      : some initialization and namelist read 
    1515   !!   ice_thd_pnd           : main calling routine 
    1616   !!---------------------------------------------------------------------- 
    17    USE phycst           ! physical constants 
    18    USE dom_oce          ! ocean space and time domain 
    19    USE ice              ! LIM-3 variables 
     17   USE phycst         ! physical constants 
     18   USE dom_oce        ! ocean space and time domain 
     19   USE ice            ! sea-ice: variables 
     20   USE ice1D          ! sea-ice: thermodynamics variables 
     21   USE icetab         ! sea-ice: 1D <==> 2D transformation 
    2022   ! 
    21    USE lbclnk           ! lateral boundary conditions - MPP exchanges 
    22    USE lib_mpp          ! MPP library 
    23    USE in_out_manager   ! I/O manager 
    24    USE lib_fortran      ! glob_sum 
    25    USE timing           ! Timing 
     23   USE in_out_manager ! I/O manager 
     24   USE lib_mpp        ! MPP library 
     25   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
     26   USE timing         ! Timing 
    2627 
    2728   IMPLICIT NONE 
     
    4041#  include "vectopt_loop_substitute.h90" 
    4142   !!---------------------------------------------------------------------- 
    42    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    43    !! $Id: limdyn.F90 6994 2016-10-05 13:07:10Z clem $ 
     43   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
     44   !! $Id: icethd_pnd.F90 8420 2017-10-05 13:07:10Z clem $ 
    4445   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4546   !!---------------------------------------------------------------------- 
    4647CONTAINS 
    4748 
    48    SUBROUTINE ice_thd_pnd( kt ) 
     49   SUBROUTINE ice_thd_pnd 
    4950      !!------------------------------------------------------------------- 
    5051      !!               ***  ROUTINE ice_thd_pnd   *** 
     
    5253      !! ** Purpose :   change melt pond fraction 
    5354      !!                 
    54       !! ** Method  :   brutal force 
     55      !! ** Method  :   brut force 
    5556      !! 
    5657      !! ** Action  : -  
    5758      !!              -  
    58       !!------------------------------------------------------------------- 
    59       INTEGER, INTENT(in) ::   kt    ! number of iteration 
    60       INTEGER ::   ji, jj, jl     ! dummy loop indices 
    6159      !!------------------------------------------------------------------- 
    6260 
     
    9694      !!     
    9795      !!------------------------------------------------------------------- 
    98        
    99       WHERE ( a_i(:,:,:) > 0._wp .AND. t_su(:,:,:) >= rt0 )  
    100          a_ip_frac = rn_apnd 
    101          h_ip      = rn_hpnd     
    102          v_ip      = a_ip_frac * a_i * h_ip  
    103          a_ip      = a_ip_frac * a_i 
    104       ELSEWHERE 
    105          a_ip      = 0._wp 
    106          h_ip      = 0._wp 
    107          v_ip      = 0._wp 
    108          a_ip_frac = 0._wp 
    109       END WHERE 
     96      INTEGER  ::   ji        ! loop indices 
     97      !!------------------------------------------------------------------- 
     98      DO ji = 1, npti 
     99          
     100         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 
     101            a_ip_frac_1d(ji) = rn_apnd 
     102            h_ip_1d(ji)      = rn_hpnd     
     103            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji) 
     104         ELSE 
     105            a_ip_frac_1d(ji) = 0._wp 
     106            h_ip_1d(ji)      = 0._wp     
     107            a_ip_1d(ji)      = 0._wp 
     108         ENDIF 
     109          
     110      END DO 
    110111       
    111112   END SUBROUTINE pnd_CST 
     
    117118      !! ** Purpose    : Compute melt pond evolution 
    118119      !! 
    119       !! ** Method     : Empirical method. A fraction of meltwater is accumulated  
    120       !!                 in pond volume. It is then released exponentially when 
    121       !!                 surface is freezing. 
     120      !! ** Method     : Empirical method. A fraction of meltwater is accumulated in ponds  
     121      !!                 and sent to ocean when surface is freezing 
     122      !! 
     123      !!                 pond growth:      Vp = Vp + dVmelt 
     124      !!                    with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i 
     125      !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) 
     126      !!                    with Tp = -2degC 
    122127      !! 
    123128      !! ** Tunable parameters : (no real expertise yet, ideas?) 
     
    131136      !!     
    132137      !!------------------------------------------------------------------- 
    133        
    134       INTEGER, DIMENSION(jpij)         ::   indxi             ! compressed indices for cells with ice melting 
    135       INTEGER, DIMENSION(jpij)         ::   indxj             ! 
    136  
    137       REAL(wp), DIMENSION(jpi,jpj)     ::   zwfx_mlw          ! available meltwater for melt ponding 
    138       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zrfrac            ! fraction of available meltwater retained for melt ponding 
    139  
    140138      REAL(wp), PARAMETER ::   zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding 
    141139      REAL(wp), PARAMETER ::   zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            '' 
    142       REAL(wp), PARAMETER ::   zrexp  = 0.01_wp  ! rate constant to refreeze melt ponds 
    143140      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp ! pond aspect ratio 
    144  
    145       REAL(wp) ::   zhi               ! dummy ice thickness 
    146       REAL(wp) ::   zhs               ! dummy snow depth 
    147       REAL(wp) ::   zTp, z1_Tp        ! reference temperature 
    148       REAL(wp) ::   zdTs              ! dummy temperature difference 
     141      REAL(wp), PARAMETER ::   zTp = -2._wp !  
     142 
     143      REAL(wp) ::   zfr_mlt            ! fraction of available meltwater retained for melt ponding 
     144      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding 
     145      REAL(wp) ::   z1_Tp        ! reference temperature 
    149146      REAL(wp) ::   z1_rhofw          ! inverse freshwater density 
    150147      REAL(wp) ::   z1_zpnd_aspect    ! inverse pond aspect ratio 
    151       REAL(wp) ::   zvpold            ! dummy pond volume 
    152  
    153       INTEGER  ::   ji, jj, jl, ij    ! loop indices 
    154       INTEGER  ::   icells            ! size of dummy array 
     148      REAL(wp) ::   zvpnd            ! dummy pond volume 
     149      REAL(wp) ::   zfac, zdum 
     150 
     151      INTEGER  ::   ji        ! loop indices 
    155152      !!------------------------------------------------------------------- 
    156153      z1_rhofw       = 1. / rhofw  
    157154      z1_zpnd_aspect = 1. / zpnd_aspect 
    158       zTp            = -2.  
    159155      z1_Tp          = 1._wp / zTp  
    160156 
    161       a_ip_frac(:,:,:) = 0._wp 
    162       h_ip     (:,:,:) = 0._wp 
    163  
    164       !------------------------------------------------------------------ 
    165       ! Available melt water for melt ponding and corresponding fraction 
    166       !------------------------------------------------------------------ 
    167  
    168       zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp )        ! available meltwater for melt ponding 
    169  
    170       ! NB: zwfx_mlw can be slightly negative for very small values (why?) 
    171       ! This can in some occasions give negative 
    172       ! v_ip in the first category, which then gives crazy pond 
    173       ! fractions and crashes the code as soon as the melt-pond 
    174       ! radiative coupling is activated 
    175       ! if we understand and remove why wfx_sum or wfx_snw could be 
    176       ! negative, then, we can remove the MAX 
    177       ! NB: I now changed to wfx_snw_sum, this may fix the problem.  
    178       ! We should check 
    179  
    180       zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:)   
    181  
    182       DO jl = 1, jpl    
    183  
    184          !------------------------------------------------------------------------------ 
    185          ! Identify grid cells where ponds should be updated (can probably be improved) 
    186          !------------------------------------------------------------------------------ 
    187          indxi(:) = 0 
    188          indxj(:) = 0 
    189          icells   = 0 
    190          DO jj = 1, jpj 
    191             DO ji = 1, jpi 
    192                IF ( a_i(ji,jj,jl) > epsi10 ) THEN 
    193                   icells = icells + 1 
    194                   indxi(icells) = ji 
    195                   indxj(icells) = jj 
    196                ENDIF 
    197             END DO 
    198          END DO 
    199  
    200          DO ij = 1, icells 
    201  
    202             ji = indxi(ij) 
    203             jj = indxj(ij) 
    204  
    205             zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    206             zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
    207  
    208             IF ( zhi < rn_himin) THEN   !--- Remove ponds on thin ice if ice is too thin 
    209  
    210                a_ip(ji,jj,jl)      = 0._wp                               !--- Dump ponds 
    211                v_ip(ji,jj,jl)      = 0._wp 
    212                a_ip_frac(ji,jj,jl) = 0._wp 
    213                h_ip(ji,jj,jl)      = 0._wp 
    214  
    215                IF( ln_pnd_fwb )   wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + v_ip(ji,jj,jl)  !--- Give freshwater to the ocean 
    216  
    217             ELSE                        !--- Update pond characteristics 
    218  
    219                !--- Add retained melt water to melt ponds 
    220                ! v_ip should never be positive, otherwise code crashes 
    221                ! MV: as far as I saw, UM5 can create very small negative v_ip values 
    222                ! hence I added the max, which was not required with Prather (1 yr run) 
    223                v_ip(ji,jj,jl) = MAX( v_ip(ji,jj,jl), 0._wp ) + zrfrac(ji,jj,jl) * z1_rhofw * zwfx_mlw(ji,jj) * a_i(ji,jj,jl) * rdt_ice 
    224  
    225                !--- Shrink pond due to refreezing 
    226                zdTs           = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. ) 
    227  
    228                zvpold         = v_ip(ji,jj,jl) 
    229  
    230                v_ip(ji,jj,jl) = v_ip(ji,jj,jl) * EXP( zrexp * zdTs * z1_Tp ) 
    231  
    232                !--- Dump meltwater due to refreezing ( of course this is wrong 
    233                !--- but this parameterization is too simple ) 
    234                IF ( ln_pnd_fwb ) THEN 
    235                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice 
    236                ENDIF 
    237                 
    238                a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) ) 
    239                !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative 
    240  
    241                h_ip(ji,jj,jl)      = zpnd_aspect * a_ip_frac(ji,jj,jl) 
    242  
    243                a_ip(ji,jj,jl)      = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl) 
    244  
     157      DO ji = 1, npti 
     158         !                                                        !--------------------------------! 
     159         IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  ! 
     160            !                                                     !--------------------------------! 
     161            !--- Remove ponds on thin ice and send water into the ocean 
     162            IF( ln_pnd_fwb ) THEN 
     163               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + h_ip_1d(ji) * a_ip_1d(ji) * rhofw * r1_rdtice 
    245164            ENDIF 
    246  
    247          END DO 
    248  
     165            a_ip_1d(ji)      = 0._wp 
     166            a_ip_frac_1d(ji) = 0._wp 
     167            h_ip_1d(ji)      = 0._wp 
     168            !                                                     !--------------------------------! 
     169         ELSE                                                     ! Case ice thickness >= rn_himin ! 
     170            !                                                     !--------------------------------! 
     171            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step 
     172            zvpnd       = v_ip_1d(ji) 
     173 
     174            ! available meltwater for melt ponding [m, >0] and fraction 
     175            zdv_mlt = -( dh_i_surf(ji)*rhoic + dh_s_mlt(ji)*rhosn ) * z1_rhofw * a_i_1d(ji) 
     176            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  !!clem CICE doc 
     177!!!            zfr_mlt = zrmin + zrmax * a_i_1d(ji)  !!clem Holland paper  
     178 
     179            !--- Pond gowth 
     180            ! v_ip should never be negative, otherwise code crashes 
     181            ! MV: as far as I saw, UM5 can create very small negative v_ip values (not Prather) 
     182            v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt ) 
     183 
     184            IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN 
     185               ! melt ponds mass flux (<0) 
     186               zfac = zfr_mlt * zdv_mlt * rhofw * r1_rdtice 
     187               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
     188 
     189               ! adjust ice/snow melting (>0) to take into account ponding 
     190               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 
     191               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
     192               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
     193            ENDIF 
     194 
     195            !--- Pond contraction (due to refreezing) 
     196            v_ip_1d(ji)      = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 
     197 
     198            ! --- Set new pond area and depth assuming linear relation between h_ip and a_ip_frac 
     199            !        h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i 
     200            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) ) 
     201            !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative 
     202            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 
     203            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji) 
     204 
     205            !! clem: there is no clever way to conserve mass here... 
     206             
     207!            IF( ln_pnd_fwb ) THEN 
     208!               ! melt ponds mass flux (>0 when ponds shrink) 
     209!               zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw * r1_rdtice 
     210!               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac!! 
     211! 
     212!               ! adjust ice/snow 
     213!               zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw / ( a_i_1d(ji)*h_s_1d(ji)*rhosn + a_i_1d(ji)*h_i_1d(ji)*rhoic )! 
     214! 
     215!               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + h_s_1d(ji) * zfac * a_i_1d(ji) * rhosn * r1_rdtice 
     216!               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     + h_i_1d(ji) * zfac * a_i_1d(ji) * rhoic * r1_rdtice 
     217!                
     218!               !!h_s_1d(ji) = h_s_1d(ji) * ( 1._wp + zfac ) 
     219!               !!h_i_1d(ji) = h_i_1d(ji) * ( 1._wp + zfac )! 
     220! 
     221!            ENDIF 
     222 
     223         ENDIF 
    249224      END DO 
    250  
    251       !--- Remove retained meltwater from surface fluxes  
    252  
    253       IF ( ln_pnd_fwb ) THEN 
    254  
    255          wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) )  
    256          wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) ) 
    257  
    258       ENDIF 
    259225 
    260226   END SUBROUTINE pnd_H12 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90

    r8597 r8623  
    164164            ! mass flux from ice/ocean 
    165165            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
    166                &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj)  
    167  
    168             IF ( ln_pnd_fwb )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj) 
     166               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 
    169167 
    170168            ! add the snow melt water to snow mass flux to the ocean 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

    r8597 r8623  
    255255      !!------------------------------------------------------------------- 
    256256      ! 
    257       v_i (:,:,:) = h_i(:,:,:) * a_i(:,:,:) 
    258       v_s (:,:,:) = h_s(:,:,:) * a_i(:,:,:) 
    259       sv_i(:,:,:) = s_i(:,:,:) * v_i(:,:,:) 
     257      v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:) 
     258      v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:) 
     259      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:) 
     260      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
    260261      ! 
    261262   END SUBROUTINE ice_var_eqv2glo 
     
    512513      ! open water = 1 if at_i=0 
    513514      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp 
     515 
    514516      ! 
    515517   END SUBROUTINE ice_var_zapsmall 
Note: See TracChangeset for help on using the changeset viewer.